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Chapter 1 

Introduction 



Gravitational wavea-^J one of the more discussed predictions of Einstein's General The- 
ory of Relativity may be detected within the next few years. General Relativity predicts 
that massive accelerating and rotating astrophysical objects emit gravitational waves 
which propagate through space with the speed of light. Gravitational waves are most 
simply thought of as ripples in the space-time fabric, their effect being to change the 
separation of adjacent masses on earth or in space; this tidal effect is the basis of all 
present detectors. Sources such as interacting black holes, coalescing compact binary 
systems, supernova explosions and pulsars are all possible candidates for detection; 
observing signals from them will significantly boost our understanding of the Universe. 

The first gravitational wave detectors were based on the effect of these tidal forces 
on the fundamental resonant mode of aluminum bars at room temperature. Initial 
instruments were constructed by Joseph Weber [3]. Following the lack of confirmed 
detection of signals, aluminum bar systems operated at and below the temperature of 
liquid helium were developed and work in this area is still underway J5J El |7]. How- 
ever the most promising design of gravitational wave detectors, offering the possibility 
of very high sensitivities over a wide range of frequency, uses widely separated test 
masses freely suspended as pendulums on earth [8j or in a drag free craft in space; 
laser interferometry provides a means of sensing the motion of the masses produced 
as they interact with a gravitational wave. Ground based detectors of this type have 

*In this thesis we decided to present only the subjects which are strictly related to the 
research topics of our work. The reader is addressed to references [T] [2J [3] for an introduction 
to the subject or an in-depth treatment of gravitational waves physics. 
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the sensitivity to observe sources whose radiation is emitted at frequencies above a few 
Hz, and space borne detectors will be developed for implementation at lower frequen- 
cies. Already gravitational wave detectors of long baseline are operating in a number 
of places around the world; in the USA (LIGO project [9] led by a Caltech/ MIT con- 
sortium) , in Italy (VIRGO project |10j . a joint Italian/French venture), in Germany 
(GEO 600 project [TT], a UK/German collaboration) and in Japan (TAMA 300 project 
|12j ) . LISA [13] , a space-based detector proposed by a collaboration of European and 
US research groups, is one of the most challenging large scale experiment of the next 
future. This detector array should have the capability of detecting gravitational wave 
signals from many astrophysical events in the Universe, providing unique information 
on testing aspects of General Relativity and opening up a new field of astronomy. 

Gravitational wave strengths are characterized by the gravitational wave amplitude 
h, given by h = 2AL/L, where AL is the change in separation of two masses a distance 
L apart. Unlike electromagnetism, gravitational radiation field is quadrupole in nature 
and this shows up in the pattern of the interaction of the waves with matter. The prob- 
lem for the experimental physicist is that the predicted magnitudes of the amplitudes 
or strains in space in the vicinity of the earth caused by gravitational waves even from 
the most violent astrophysical events are extremely small, of the order of 10 -21 or lower 
(for a review of the gravitational wave sources and expected signal strength see [14 ). 
Indeed current theoretical models on the event rate and strength of such events sug- 
gest that in order to detect a few events per year, from coalescing neutron star binary 
systems for example, an amplitude sensitivity close to 10~ 22 over timescales as short as 
a millisecond is required. If the Fourier transform of a likely signal is considered, it is 
found that the energy of the signal is distributed over a frequency range or bandwidth 
which is approximately equal to 1/timescale. Thus detector noise levels must have 
an amplitude spectral density lower than ~ 10~ 23 (Hz)~ 1 / 2 over the frequency range 
of the signal. The weakness of the signal means that limiting noise sources like the 
thermal motion of molecules in the detector (thermal noise) , seismic or other mechani- 
cal disturbances, and noise associated with the detector readout, whether electronic or 
optical, must be reduced to a very low level. For signals above ~ 10Hz ground based 
experiments are possible, but for lower frequencies where local fluctuating gravitational 
gradients and seismic noise on earth become a problem, it is best to consider developing 
detectors for operation underground [15 or in space [T3] . 
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Gravitational wave detector using laser interferometry, offers the possibility of very 
high sensitivities over a wide range of frequency. It uses test masses which are widely 
separated and freely suspended as pendulums to isolate against seismic noise; laser 
interferometry provides a means of sensing the motion of these masses produced as 
they interact with a gravitational wave. 

This technique is based on the Michelson interferometer and is particularly suited 
to the detection of gravitational waves as they have a quadrupole nature. Waves prop- 
agating perpendicular to the plane of the interferometer will result in one arm of the 
interferometer being increased in length while the other arm is decreased and vice versa. 
The induced change in the length of the interferometer arms results in a small change 
in the intensity of the light observed at the interferometer output. The sensitivity of 
an interferometric gravitational wave detector is limited by noise from various sources. 
Taking this frequency dependent noise floor into account, the American LIGO detectors 



have a sensitivity (shown in Fig. 1.1), which would allow a reasonable probability for 
detecting gravitational wave sources. 

In order to observe a full range of sources and to initiate gravitational wave as- 
tronomy, a sensitivity or noise performance approximately ten times better in the mid- 
frequency range and several orders of magnitude better at 10 Hz, is desired. Such a 
performance is planned for a future LIGO upgrade, Advanced LIGO [16J. 

In general, for ground based detectors the most important limitations to sensitivity 
result from the effects of seismic and other ground-borne mechanical noise, thermal 
noise associated with the test masses and their suspensions and shot noise in the pho- 
tocurrent from the photodiode which detects the interference pattern. 

Most modern designs implement improved versions of a simple Michelson inter- 



ferometer (see Fig. 1.2). A simple Michelson interferometer has an antenna response 
function which has a maximum sensitivity for r r t = T gw /2, where r r t = 2L/c is the 
round-trip travel time for photons leaving and returning to the beamsplitter and T gw 
is the period of the gravitational wave signal (L the length of each arm). A simple 
calculation from this expression shows that for frequency between 10Hz and 1kHz, the 
optimal antenna length is of order 10 5 m to 10 7 m. This is much larger than would be 
feasible for an earth based detector. 

The situation can be helped greatly if a multi-pass arrangement is used in the arms 
of the interferometer as this multiplies up the apparent movement by the number of 
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Figure 1.1: Sensitivity of the LIGO detector HI during the science run S5. 

bounces the light makes in the arms. The multiple beams can either be separate as in 
an optical delay line, like the GEO 600 configuration, or may lie on top of each other as 
in a Fabry-Perot resonant cavity, used in all others interferometric gravitational waves 
detectors. A Fabry-Perot cavity consists of a partially transmitting input mirror and 
a high reflective rear mirror. If the length of the Fabry-Perot is adjusted to a multiple 
of the laser wavelength the cavity becomes resonant. The light power inside the cavity 
builds up and simulates the effect of sending the light forth and back multiple times. 
However, in this case the number of bounces is not a fixed quantity, but rather an 
averaged effective value. Optimally, the light should be stored for a time comparable 
to the characteristic timescale of the signal. Thus if signals of characteristic timescale 
1 ms are to be searched for, the number of bounces should be approximately 50 for an 
arm length of 4km. 

If the mirrors have low optical losses and if the rear mirror is a high reflector, most 
of the power incident to a Fabry-Perot arm cavity will be reflected back to the beam 
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Figure 1.2: Optical scheme of the LIGO GW interferometers. 

splitter. Ideally, the anti-symmetric port of the Michelson interferometer is set on a 
dark fringe to minimize shot noise [T7] . Then a differential length change induced by a 
gravitational wave will leave through the anti-symmetric port with the highest possible 
signal-to-noise ratio. This in turn means that most of the injected light will leave the 
interferometer through the symmetric port and be lost. By placing an additional par- 
tially transmitting mirror at input one can form yet another cavity, the power recycling 
cavity, and recycle most of the otherwise lost light. The interferometer response is then 
enhanced by the power recycling gain (the shot noise is reduced by the additional power 
build-up in the power recycling cavity) and this configuration is called Power-recycled 
Michelson with Fabry-Perot arm cavities. By adding a partially transmitting mirror to 
the anti-symmetric output port the gravitational wave signal can be made resonant [18 . 
This makes it possible to shape the interferometer response, so that its sensitivity is 
improved in a narrow frequency band around the signal resonance. In general, this 
means that the sensitivity outside the resonant frequency band will be worse. This 
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Figure 1.3: Advanced LIGO sensitivity design. 

is not a problem at lower frequencies where the interferometer is usually limited by 
seismic noise. If both power and signal recycling are implemented the configuration is 
called dual recycled. 

1.1 Topics discussed in this thesis 

Thermal noise of the test masses is expected to be a limiting factor in Advanced GW 
interferometers. As an example, Fig. |1.3| shows the expected sensitivity curve for Ad- 
vanced LIGO interferometer with the main noise contribution enlightened. Many re- 
search groups work on R&.D activities finalized to improve the thermal noise perfor- 
mance of next generation detectors. Some research lines deal with cryogenic temper- 
ature, other with improved or new materials, other with optical beam shaping and 
optimization of the mirror geometry and/or coating. Non Gaussian beams have been 
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1.1 Topics discussed in this thesis 



proposed years ago |19( [20] to reduce a particular type of thermal noise (substrate 
thermoelastic) . 

In this thesis we provide a quantitative analysis of the impact of non- Gaussian beams 
on different kinds of thermal noises^ We show that the mesa beam implementation 
could boost the Advanced LIGO sensitivity considerably: even with a rough estimation 
(without re-optimizing the detector for the introduction of mesa beams), the binary 
neutron star inspiral range increases from 175 Mpc to 225 Mpc. 

We illustrate the importance of uniform sampling of the mirror surface to reduce 
thermal noise and the limitation brought by the use of excited modes with nodes on 
the mirror surface. 

We developed the theory of mesa beam, in view of a future implementation in ad- 
vanced GW interferometers of the mesa beam idea, focusing on the analytical derivation 
of the quantities (beam width, divergence, M 2 factor, etc.), which are chosen as "ISO 
standard" [22] reference parameters for the characterization of an optical beam. 

We also analytically proved a new duality relation between optical cavities with non- 
spherical mirrors. This derivation provides a unique mapping between the eigenvalues 
and eigenvectors of two cavities whose mirrors shapes are related by a simple relation. 
This duality allows the direct application of beam property calculations performed in 
a case to geometries of the other configuration. 

The interest of the GW community in this new beam technology led us to the con- 
struction and testing of a prototype mesa beam Fabry-Perot cavity with mexican-hat 
mirror. Part of the work of this thesis was devoted to the development of new simula- 
tion programs of optical system^} These programs provided the theoretical expected 
behavior of our experiment, in particular cavity's modes structure and misalignments 
sensitivity to be confronted with the experimental results. We developed new simula- 
tion packages to analyze the performance of our cavity prototype with real imperfect 
mirrors, using the measured mirrors maps. The model developed can include uniform 
and non-uniform scattering and absorption losses, as well as the effects of mirror heat- 
ing A particular attention has been devoted to keep theses simulation programs very 

-'In particular we show that the thermal noise of the mirror's dielectric coating is greatly 
reduced by using flat profile beams. This resolved a question raised [3T] on the effectiveness of 
this type of beam on coating noise 

2 Available at http://www.ligo.caltech.edu/~jagresti/ 
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1.1 Topics discussed in this thesis 



easy to use and very easy to change/upgaddj by anyone involved in optics researcl 




The good agreement between theory and experiment validated the mathematical tools 
here developed thus allowing safe extrapolation to the larger optical systems needed in 
GW observatories. 

We also explored another complementary way of reducing the mirror thermal noise, 
beside the beam shaping, that is the multi-layered coating thickness optimization. We 
show it to be effective in reducing the coating noise and explore the possible implications 
for GW interferometers in terms of sensitivity. 

During this analysis we developed an independent model for the coating effective 
elastic parameters, which is based on the well understood subject of homogenization 
theory. 



1 Most of the time we used Mathematica® for its capability of analytical manipulation and 
numerical analysis or Matlab® for its ability in manipulating large matrices. 

2 The great simplicity and versatility of these programs turned to be fundamental properties 
for the usage by other RSzD research groups. 
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Chapter 2 

Paraxial beam and optical 
cavities 

2.1 Introduction 

This chapter is a brief description of the techniques used in this thesis to study parax- 
ial laser beam propagation and optical resonators. It gives an overview of the main 
concepts related to standard spherical mirrors optics, and some less known concepts 
related to generic paraxial beam and resonators which is one of the main subject of 
this work. The reader familiar with these techniques can jump to Chapter [3] 

2.2 Paraxial laser beams 

Radiation from lasers is different from conventional optical light because it is very 
close to be monochromatic. Although each laser has its own fine spectral distribution 
and noise properties, the electric and magnetic fields from lasers are considered to 
have minimal phase and amplitude variations in the first-order approximation. Like 
microwaves from a maser, electromagnetic radiation with a precise phase and amplitude 
is described most accurately by Maxwell's wave equations. For devices with structures 
that have dimensions very much larger than the wavelength, e.g. in a multimode fiber 
or in an optical system consisting of lenses, prisms or mirrors, the rigorous analysis of 
Maxwell's vector wave equations becomes very complex and tedious: there are too many 
modes in such a large space. It is difficult to solve Maxwell's vector wave equations 
for such cases, even with large computers. Even if we find the solution, it would 



9 



2.2 Paraxial laser beams 



contain fine features (such as the fringe fields near the lens) which are often of little or 
no significance to practical applications. In these cases we look for a simple analysis 
which can give us just the main features (i.e. the amplitude and phase) of the dominant 
component of the electromagnetic field in directions close to the direction of propagation 
and at distances reasonably far away from the aperture. When one deals with laser 
radiation fields which have slow transverse variations and which interact with devices 
that have overall dimensions much larger than the optical wavelength A, the fields can 
often be approximated as transverse electric and magnetic (TEM) waves. In TEM 
waves both the dominant electric field and the dominant magnetic field polarization 
lie approximately in the plane perpendicular to the direction of propagation. For such 
waves, we usually need only to solve the scalar wave equations to obtain the amplitude 
and the phase of the dominant electric field along its local polarization direction. The 
dominant magnetic field can be calculated directly from the dominant electric field. 
Under these assumption, the Maxwell equations for the beam propagation in free space 
(or homogeneous and isotropic medium), can be replaced by the scalar wave equation 
(Helmholtz equation) 

[V 2 + k 2 ]E(x, y, z) = (2.1) 

where E(x, y, z) is the phasor amplitude of a field component that is sinusoidal in 
time, V 2 is the Laplacian operator and k is the laser wavenumber. Then, since laser 
beams are usually sufficiently collimated, we can describe their diffraction properties 
using the paraxial wave approximation. If the field is expected to propagate mainly 
in the z direction, with a slow variation of amplitude and phase along the transverse 
direction, it is convenient to write the field in the following way 

E(x,y,z) = u{x,y,z)e- ikz (2.2) 

where u is the complex scalar wave amplitude, called envelope function, which de- 
scribes the transverse profile of the beam. Inserting this into the (Helmholtz) equation 



(2.1 ) we find the reduced equation 

d 2 2 <)~ . ., 
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2.3 Integral approach 



If the z dependence of the envelope function is slow compared to the optical wave- 
length and to the transverse variations of the field, we can drop the second order partial 



derivative in z in (2.3) and obtain the paraxial wave equation 



d 2 u d 2 u 2'k^ U (2 4) 

dx 2 dy 2 dz 
This equation can be written for a generic set of transverse coordinates as 

VMs,z)-M^M = (2.5) 

where s refers to the specific coordinate system, orthogonal to the z direction, and 
V 2 means the Laplacian operator in these coordinates. 



2.3 Integral approach 

Another equally valid and effective way of analyzing paraxial wave propagation is to 
employ the Huygens-Fresnel principle in the paraxial approximation. Consider an aper- 
ture at a plane zq illuminated with a light field distribution Eq(xq, yo, zq) within the 
aperture. Then for a point lying somewhere after the aperture, say at P with coordi- 
nates (x,y,z), the net field is given by adding together spherical waves emitted from 
each point Pq in the aperture. Each spherical wavelet takes on the strength and phase 
of the field at the point where it originates. Mathematically, this summation takes the 
form 



i I , ,e~ tkp 



E(x,y,z) = - E (xo,yo,zo) cos 9 dS (2.6) 

A J s Q P 

P = V( z - z o) 2 + x ) 2 + {y - yo) 2 

where 9 is the angle under which the aperture element centered at Po is seen from 
the observation point P and p is the distance between these two points. This equation 
can be rigorously derived from the Rayleigh-Sommerfeld scalar diffraction theory as 
shown in [23] and rely on the approximation p S> A. The paraxial and the Fresnel 
approximations to diffraction theory consist in approximating the obliquity factor cos 9 
by unity and, once expanded the distance p in a power series in the form 

, (x - x ) 2 + (y- y ) 2 
p = z-z + — r + . . . (2.7 

2(Z- Zq) 
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2.3 Integral approach 



replacing the denominator of (2.6) by simply z — zq and retaining the quadratic 



terms in the exponent of the phase shift factor e lkp . With these approximations we 



obtain the paraxial approximation of the Huygens-Fresnel integral (2.6) 



E(x,y,z) = — - E (x ,y ,z )e 2 ^~*o) dx dy (2.8) 

A{Z — Zq) J 

The integral is extended to the entire plane zq with the assumption that the input 
field distribution Eq{xq,uq, zq) is a paraxial optical beam which transverse profile going 



to zero outside the region Sq. The integral (2.8) provides a way to propagate an 
arbitrary optical wavefront from an input plane zq to any later plane z. In the following 
we will assume zq = 0. It is useful to note that this integral has exactly the form of a 
convolution product in the x, y coordinates between the input field Eq(xo, yo, zq = 0) 
and the paraxial diffraction kernel h(x, y, z) 

A s~i ikz 2 i 2 

h(x,y,z) = —— e 2, (2.9) 

AZ 



Using the 2D Fourier transform the integral equation ( 2.8 ) can be transformed into a 
simple algebraic product. Consider a complex function (e.g. an optical wave amplitude) 
f(x, y) of two real variables, of integrable square modulus. Its two dimensional Fourier 
transform is defined by 



f(p,q)= I dxdye-^ x+ ^f(x,y) (2.10) 
and the reciprocal transform by 



R 



2 



f(x,y) = ^J 2 dpdqe^ x +^f(p,q) (2.11) 
The Fourier transform of the propagator h can be easily computed 

h(p,q,z) = e - lkz e izE ^ (2.12) 
and the propagation equation for the transformed field become 

E(p, q, z) = h{p, q, z)Eo{p, q, 0) (2.13) 
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2.3 Integral approach 



This equation can be retrieved starting from the Helmholtz equation in the paraxial 



approximation. Taking the Fourier transform of the equation (2.4) with respect to x,y 



we obtain the following equation for the Fourier transformed envelope function 

[2ikd z + {p 2 +q 2 )]u{p,q,z) =0 (2.14) 
the solution of which is of the form 

u(p, q, z) = u(p, q, 0)e tzIL ^k~ (2.15) 

in which we recover the propagator ( 2.12| ) once we reintroduce the phase factor 
e -ikz £ Qr |- ne com pi e t e field. 

2.3.1 Paraxial plane waves decomposition 



The scalar wave equation (2.1 ) gives us a formal method for propagating an optical wave 
forward in space. The main concept in Fourier optics is that we use a Fourier transform 
to break up an arbitrary field into plane- wave components. Then we propagate the 
plane waves components, and the linearity of the propagation equations allows the 
reconstruction of the propagated field. Let's see how this works. Recall the scalar 
plane wave solution to the wave equation: 

E(r) = E e~ ik - r = E Q e- i{k * x+k y y+k * z) (2.16) 

where the wave vector k = (k x ,k y ,k z ) indicates the propagation direction of the 
wave. Lets let the optical axis lie along the z-axis. Then the wave vector k has a 
direction determined by its angles with respect to the z-axis, 

kx \ i/i __• ( 



6 X = arccos yjr j an< ^ ®v = arcsni J (2-17) 
where k x and k y are the transverse spatial frequencies. 

That is, in the z = plane, the electric field has harmonic spatial dependence in 
the x- and y-directions: 

E(x, y,z = 0) = E e- i{kxX+k y y) (2.18) 
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2.3 Integral approach 



This is the central point of Fourier optics: harmonic phase variation in the (x, y) 
plane (say, at z = 0) corresponds to a plane wave in a particular direction determined 



by (2.17). That is, the frequency of the harmonic variation determines the direction 
of propagation. A general field profile E(x, y) at z = can thus be written as a 
superposition of plane waves via the Fourier transform: 



E(x,y) 



1 

4^2 JR2 

where the spatial frequency distribution is given by 



dk x dk y e i{k * x+k y y) E(k x ,k y ) 



(2.19) 



E(k X i ky 



R 2 



dxdye- i{k * x+k y y) E{x,y) 



(2.20) 



Having transformed E(x,y) into the spatial-frequency domain, we can propagate 
this spatial-frequency distribution forward to any other plane z by multiplying it by 
the phase shift factor e~ lkzZ , where k z = . k 2 — k x — ky 



E(k x , k y , z) = E(k x , k y )e 



-ik z z 



(2.21) 



However in the paraxial approximation, since k x<u <C k , the longitudinal component 
k 7 can be written in the form 



I/-- _ . _ £ r"*-^ l/i 

hi fVg, rvy ~ rv 



1 / kx H~ k, 



(2.22) 



Thus, the free-space transfer function becomes 



-ik z z 



(T lkz e % 



(2.23) 



where the first factor is an overall phase factor corresponding to plane-wave propa- 
gation along the optical axis, and the second factor generates the evolution of the spatial 



profile. Combining Eqs. (2.13) and (2.23) we obtain the expression for the propagation 



of the spatial-frequency components of the paraxial beam which is equal to Eq. (2.13), 
obtained by applying Fourier theorems to the Huygens-Fresnel integral. 
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2.4 Gaussian beam 



2.4 Gaussian beam 

The Gaussian beam is the simplest model of a directed beam that satisfies the Helmholtz 



equation in the paraxial approximation (2.4). It also turns out that the outputs of 



spherical mirror resonators and lasers are often Gaussian beams (approximatively), 
and therefore the theory of Gaussian beams is widely used in gravitational waves in- 
terferometric detectors. In most applications the standard Gaussian beam is expressed 
as 



E G (r) = E c 



WQ 

w(z) 



exp 



w 2 (z) 



exp 



-ikz + i arctan 



zo 



exp 



-ik 



2R{z) 



(2.24) 



In this expression, we have used polar coordinates (r, z) with r = \J x 2 + y 2 ; Eq 
is an overall field-amplitude constant; zq is a constant called the Rayleigh length (or 
"Rayleigh range"); 



w 



Xzo 

7T 



(2.25) 



is the beam waist parameter or minimum beam radius 



w(z) = w Q \ 1 + 



is the beam radius as a function of z, and 



zo 



(2.26) 



R{z) 



1 + 



50 y 

z J 



(2.27) 



is the radius of curvature of the wave front. The intensity is the square modulus of 



Eq. (2.24) 



I(r,z) = \E Q f 



w 
w(z) 



exp 



2r 2 



w 2 {z) 



(2.28) 



Clearly, the intensity falls off in the radial direction like a Gaussian function, hence 
the name Gaussian beam. Also, note that the beam spot size w{z) correspond to the 
distance from the optic axis to the point of 1/e 2 attenuation of the beam intensity (1/e 
in electric field amplitude). 
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2.4 Gaussian beam 



The beam radius w{z) traces out a hyperbolic curve in z. Near the focus at z = 0, 
the spot size achieves its minimum value wq- At large distances from the focus, the 
hyperbola approaches its asymptotes, given by {wq/zq)z. The Rayleigh length zq marks 
the crossover between these two regimes; thus, 



w(z) ~ wo for \z\ <§; Zq 
w(z) ~ ( — ) z for \z\ 3> Zq 

for the regions near and far away from the focus, respectively. From this we see 
that in the far field, the beam propagates in the form of a cone of half angle 9o, called 
beam divergence, given by 

#o~tan# = — = — (2.29) 

Z TTW 

where we have used the small angle approximation for the paraxial Gaussian beam. The 
region in which a gaussian beam can be considered collimated is roughly 2zq around 
the focal plane z = 0. Recall that the longitudinal phase factor has the form 



exp 



-ikz + i arctan 



(2.30) 



z 

. z o. 

The first term in the phase is simply the phase of a plane wave ikz propagating in the 
same direction and with the same optical frequency as the Gaussian beam. The second 
term is called the Gouy phase shift and represents a small departure from planarity. 
The longitudinal phase is dominated by the ikz term, but the Gouy term is important 
as well. It represents a phase retardation compared to the plane wave. Because of the 
arctan form, the retardation amounts to a total of tt in phase over all z. Gouy effects 
are generic to focusing-beam-type solutions to the wave equation. As we will see, the 
Gouy phase is important in computing the resonant frequencies of optical resonators. 



Examining the expression of the radius of curvature of the wavefront (2.27) we 
can see that at the waist z = 0, the wavefront is flat, and that in the far field region 
(\z\ 3> zo) the radius increases as R(z) « z, i.e., the gaussian beam becomes essentially 
like a spherical wave centered at the beam waist. The minimum radius of curvature 
occurs for the wavefront at a distance from the waist given by z = Zq, with a radius 
value R = 2zq. In a cavity, the boundary conditions imposed by the cavity mirrors 
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2.5 ABCD transformation 



require that the curvature of the spherical mirrors and the curvature of the wave fronts 
match. This allows the wave to map back on itself. This is one reason why Gaussian 
beams can exist in resonators as we will see later. 

Despite the complex form of the Gaussian beam, relatively little information is 
needed to completely specify it. For example, if we know where z = is located, the 
value of wo, and the optical wavelength A all the other parameters are uniquely fixed. 
Alternately, it is sufficient to know u>o and R{z) at some distance z, or it is sufficient 
to know w(z) and R{z) at some distance z. It is useful the definition of a complex 
parameter for the gaussian beam, called g-parameter 

q(z)=z + iz (2.31) 



Using the definition (2.27) and (2.26) it is easy to calculate 1/q in the more useful 
form 

1 1 A , 

(2.32) 



q(z) R{z) ttw 2 {z) 
The propagation of a gaussian beam in free space can be simply computed using 
the q-parameters; the propagation law between two planes along the optical axis z for 
q(z) is established as 



q(z 2 ) = q(zi) + Z2- z 1 (2.33) 

This result is a particular case of the propagation low of a gaussian beam through 
various optical structures modelled with the so-called "ABCD" formalism of geometric 
optics which we will describe briefly in the next section. 



2.5 ABCD transformation 

Matrix optics has been well established a long time ago. Within the paraxial approach, 
it provides a modular transformation describing the effect of an optical system as 
the cascaded operation of its components. Then each simple optical system is given 
by its matrix representation. Before presenting the results of the application of the 
matrix optics to the Gaussian beam transformation, we need to analyze the basis of 
this approach (e.g., see Chapter 15 of Ref. [24 J. In paraxial optics, the light is presented 
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2.5 ABCD transformation 



as ray trajectories that are described, at a given meridional plane, by its height and 
its angle with respect to the optical axis of the system. These two parameters can be 
arranged as a column vector. The simplest mathematical object relating two vectors 
(besides a multiplication by a scalar quantity) is a matrix. In this case, the matrix is a 
2x2 matrix that is usually called the ABCD matrix because its elements are labelled 
as A, B, C, and D. If the input plane and the output plane are in the same optical 
medium, then the determinant of the ABCD matrix is unity, AD — BC = 1, and the 
matrix is called unimodular. In the absence of loss, the matrix elements are real, they 
are complex otherwise. 

The relation can be written as: 



x 2 



A B 
C D 



(2.34) 



where the column vector with subindex 1 stands for the input ray, and the subindex 
2 stands for the output ray. An interesting result of this previous equation is obtained 
when a new magnitude is defined as the ratio between height and angle. 

Optkal 
System 



A 


B 


C 


D 




Optical Axis 



Input Plane 



Output Plane 



Figure 2.1: Ray transformation trough an ABCD system 



From Fig. 2.1 this parameter coincides with the distance between the ray-optical 
axis intersection and the position of reference for the description of the ray. This 
distance is interpreted as the radius of curvature of a wavefront departing from that 
intersection point and arriving to the plane of interest where the column vector is 
described. When this radius of curvature is obtained by using the matrix relations, the 
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2.5 ABCD transformation 



following result is found: 

ARi + B . . 

fl2 = cfihn? < 2 - 35) 

This expression is known as the ABCD law for the radius of curvature. It relates the 
input and output radii of curvature of a spherical wavefront passing through an optical 
system described by its ABCD matrix. In recent years, a very useful generalized form 
of paraxial beam propagation has been developed, in which a generalized Huygens type 
of integral, describes paraxial wave propagation trough cascade sequences of optical 



elements, in terms of ABCD matrices. Consider an optical system (Fig. 2.1 ) described 
by an ABCD matrix with real elements and illuminated by a beam with a transverse 
amplitude distribution E\{x\,y{) at the input plane. The amplitude distribution at 
the output plane can be written as 

E 2 (x 2 ,y 2 ) = ^r-^ f J Bi(xi,yi)e- i ^[ A( ^ + ^ ) ~ 2(xiX2+yi?/2)+D(:i -2 + ^)]dx 1 dy 1 (2.36) 
XB J 

or for an optical system with cylindrical symmetry 



E 2 (r 2 ,> 



Hi 



-ikL 



XB 



' m [ Ar i - 2ri ra cos (0i _ ^ 2 ) + Dr 2 



r\dr\d<f>\ 



(2.37) 



If we now consider the transformation of a gaussian beam through a paraxial system 
described by an ABCD matrix, it is possible to prove (ref..) that the complex q- 



parameter introduced in (2.32) transforms according to the simple relation 



<?2 



A qi + B 



(2.38) 



Cq x + D 

where q± and q 2 are the ^-parameters before and after the optical system, respec- 



tively. The equation (2.38) has the same form of the equation (2.35) for the transforma- 



tion of the radius of curvature of a spherical wave passing through the optical elements; 
the q-parameter is also called complex radius of curvature for a gaussian beam. The 
results of the application of the ABCD law can be written in terms of the real radius 
of curvature R and the Gaussian width w by properly taking the real and imaginary 
parts of the resulting complex radius of curvature. When a Gaussian beam propagates 
along an ABCD optical system, its complex radius of curvature changes according to 
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2.6 Higher order modes 



the ABCD law. The new parameters of the beam are obtained from the value of the 
new complex radius of curvature. However, there exists an invariant parameter that 
remains the same throughout ABCD optical systems. This invariant parameter is the 
product of the minimum beam width (wq) with the divergence of the beam 

w 6 = A (2-39) 

7T 

Using this relation, we can conclude that a good collimation (very low value of the 
divergence) will be obtained when the beam is wide. On the contrary, a high focused 
beam will be obtained by allowing a large divergence angle. 



2.6 Higher order modes 



In section 2.4 only one solution of (2.5), was discussed, i.e., a light beam with the 
property that its intensity profile in every beam cross section is given by the same 
function, namely, a Gaussian. The width of this Gaussian distribution changes as the 



beam propagates along its axis. There are other solution of (2.5 ) with similar properties, 
and they are discussed in this section. These solutions form a complete and orthogonal 
set of function and are called the modes of propagation. Every arbitrary distribution 
of monochromatic light can be expanded in terms of these modes. This procedure is 
largely used in the modal analysis of the light circulating in the gravitational waves 
interferometers. 

a) Modes in cartesian coordinates 

Using cartesian coordinates, a more general solution of the paraxial wave equation 
is given by the Hermite-Gaussian (HG) beam which can be written as 



HG mn {x,y,z) 



7r2 m+n m\n\w 2 (z) 



exp 



-ikz + i(m + n + 1) arctan 






x 2 + y 2 ' 


I exp 


w 2 (z) 





-ik 



x 2 + y 2 
2R{z) 



(2.40) 



where the functions H m {X) are the Hermite polynomials of order m and the pa- 
rameters w(z),R(z) and zq are the same as for the lowest-order gaussian mode as given 



in Sec. 2.4 In general, the intensity pattern of the HG myn mode, Fig. 2.2 has m dark 
bands across the x-direction and n dark bands across the y-direction, corresponding to 
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2.6 Higher order modes 



the zeros of the Hermite polynomials. Alternately, the intensity pattern has a grid of 
m + 1 bright spots in the x-direction and n + 1 in the y-direction. 



HGnn HGni HGn 






HGio HGn HG12 

4^ 4^ 4^ 

HG20 HG21 HG22 

Figure 2.2: Hermite Gauss modes 
b) Modes in cylindrical coordinates 

An alternative but equally valid family of solution to the paraxial wave equation 
( |2.5[ ) can be written in cylindrical coordinates (r, <ft, z) and are called Laguerre Gaussian 
(LG) beams 



LG prn (r, <f>, 



1 4pl ( V2r 

7r(l + 6 m0 )(m + p)\ \w(z) 



t m 

V 



2r 2 



w 2 (z) 



exp 



w 2 (z) 



w(z) 



exp 



-ikz + i{2p + m + 1) arctan 



cos(m^>) exp 



-ik 



2R(z) 
(2.41) 



where the integer p > is the radial index and the integer m is the azimuthal mode 
index (see Fig. 



2.3 



for the Intensity distribution); the L™ are the generalized La- 
guerre polynomials and all other quantities R(z), w(z), zq are exactly the same as in 
the Hermite-gaussian case. 
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2.7 Generic paraxial beam and ABCD transformation 



LGn 



LGn 



LGn 




LGi 



LGi 



LGi 




LG, 





LG21 





LG22 



Figure 2.3: Laguerre Gauss modes 



2.7 Generic paraxial beam and ABCD transformation 



The propagation and transformation of Gaussian, Hermite-Gaussian and Laguerre- 



Gaussian light beams through paraxial systems are governed by the ABCD law (2.38). 



However, since most of the work in this thesis deals with non-Gaussian (non-HG and 
non-LG) beams, it is useful to introduce a formalism for the analysis of a non-Gaussian 
beam by using an extension of the complex beam q-parameter. We are going to intro- 
duce the most characteristic parameters for a generic paraxial beam, defined in terms 
of the moments of the intensity distribution and its Fourier transform. These parame- 
ters have been introduced in |25( |26| [27] . In the case of orthogonal astigmatic beams, 
the field is separable in two orthogonal one-dimensional components. If the amplitude 
distribution in one of these transversal directions is denoted by ip(x), then the width 
w(i/j) of the beam is given by 



where 



u#) = 2 I x (V) 

/•CO 

= \ \ifj( x )\ 2 dx 



(2.42) 
(2.43) 
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2.7 Generic paraxial beam and ABCD transformation 



is the total intensity in the transverse direction and 

1 f°° 

= Ju^J \^)\ 2 xdx (2.44) 

is the mean value of the transversal position of the beam. When x(tp) is zero, we 
will say that the beam is on axis. It is easy to check that in the case of a Gaussian 
distribution, the width is the Gaussian width defined in the previous sections. 

As we saw in the definition of the divergence for Gaussian beams, the divergence is 
related to the spreading of the beam along its propagation. This concept is described 
analytically by the Fourier transform of the amplitude distribution, i.e., also named as 
the angular spectrum. In this section the Fourier transform of the amplitude distribu- 



tion is defined in a slightly different way with respect to (2.10); introducing a 2ir factor 
in the exponent we have 

/oo 
ifj(x)e- i2 ^ x dx (2.45) 
-oo 

where £ is the transverse spatial frequency that is related to the angle by means 
of the wavelength. The angular width that is taken as the divergence of the beam is 
related to the mean-square deviation of the Fourier transform of the amplitude 

flbWO = 2A f " £ 2 (V0J (2-46) 

where due to the Parseval's theorem = I{ip) and 

~ 1 r°° ~ 

m = — =- / IVf (2-47) 

I(lp) J-oo 

is the mean value of the transversal spatial frequency which is related to the angle be- 
tween the optical axis and the direction of the propagation of the beam by a = 
The mean value of the position, x(ip), and the slope of the beam, a transform according 



to the geometrical optics rules (2.34) when the beam passes trough the ABCD system. 

Another parameter defined in analogy with the Gaussian beam case is the radius 
of curvature. For totally coherent laser beams, it is also possible to define an effective 
or generalized radius of curvature for arbitrary amplitude distributions. This radius of 
curvature is the radius of the spherical wavefront that best fits the actual wavefront of 
the beam. This fitting is made by weighting the departure from the spherical wavefront 
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2.7 Generic paraxial beam and ABCD transformation 



with the irradiance distribution. The analytical expression for this radius of curvature 
can be written as follows: 



t°° f d ^%*( x ) _ ) [ x _ xW]dx (2.48) 



R(ip) Trl{ijj)w 2 {il)) J_ 00 \ dx dx 

By using the previous definitions, it is possible to introduce a generalized complex 
radius of curvature as follows: 



1 - 1 .-./«<« 1 (2.49) 



Now the transformation of the complex radius of curvature can be carried out by apply- 



ing the ABCD law, as it can be proved by using the generalized Huygens integral (2.36) 



or (2.37), which relates the output beam to the input one and the ABCD elements. It 
is important to note that there are three parameters involved in the calculation of the 
generalized complex radius of curvature: w 2 (ij:), 0q(iP) and R(ip). The application of 
the ABCD law provides two equations: one for the real part, and one for the imaginary. 
Therefore we will need another relation involving these three parameters to solve the 
problem of the transformation of those beams by ABCD optical systems. This third 
relation is given by the invariant parameter M 2 called beam propagation /actorQ 

For the Gaussian beam case, we have found a parameter that remains invariant 
through ABCD optical systems. Now in the case of totally coherent non-Gaussian 
beams, we can define a new parameter that will have the same properties. It will be 
constant along the propagation through ABCD optical systems. Its definition in terms 
of the previous characterizing parameters is: 



M2 = l wW f 2 ^ ) -RW) (2 - 50) 
This invariance, along with the results obtained from the ABCD law applied to the 
generalized complex radius of curvature, allows to calculate the three resulting param- 
eters for an ABCD transformation. The value of the M 2 parameter has an interesting 
meaning. It is related to the divergence that would be obtained if the beam having 
an amplitude distribution if] is collimated at the plane of interest. The collimation 
1 Sometimes it is also called beam quality factor. 
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2.7 Generic paraxial beam and ABCD transformation 



should be considered as having an effective, or generalized, radius of curvature equal to 
infinity. From the definition of R{ip), this is an averaged collimation. The divergence of 
this collimated beam is the minimum obtainable for such a beam having a generalized 
width of w(ip). The propagation factor M 2 can be directly related to the product of 
the minimum width (waist size) of the beam (defined as variance in the x coordinate) 
times the divergence (defined as variance in the £ coordinate). The value of M 2 for a 



Gaussian beam is 1, which follows directly by (2.39). It is not possible to find a lower 
value of the M 2 for actual, realizable beams. This property, along with its definition 
in terms of the variance in x and £, resembles very well an uncertainty principle. Us- 
ing these definitions, the transformation of the beam width, divergence and radius of 
curvature of a generic paraxial beam by a real ABCD system from an input plane to 
an output plane (for simplicity in the following formulas we will assume that the beam 
is on axis and that the slope is zero at the input plane) can be written as 



W 2 (lp2) = W 2 (t{)l) 

1 w 2 (^) 



A + 
C + 
A + 



B 


2 


+ B 2 


M 4 X 2 






■K 2 W 2 {lj) 


D 
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+ D 2 


M 4 A 2 


i?(Vi)J 




TT 2 W 2 (lfj 


B 




c + 


D 


R(^i) 







+ BD 



M 4 A 2 



TT 2 W 2 (lpl)'W 2 (lJj2 



(2.51) 
(2.52) 
(2.53) 



In particular, if z = is the plane of the smallest width, the free propagation of the 
beam (A = D = 1,C = and B = z) produces these equations 



w 2 (ip z 
R{^z 



w 2 (ipo 



1 + 



1 + z 



2 4W) 
W 2 (lp ] 



w 2 (i/j 



1 + 



w 2 (ip ) 



?o 2 (^o) 



= z 


H?)1 







z 

ZR 



where we have introduced a generalized Rayleigh distance zr 



TTW 2 (lpo) 

M 2 \ ■ 



(2.54) 
(2.55) 
Since in 



the next chapter we will deal with cylindrical symmetric laser beam, it is useful to write 
the definitions of width, divergence, and radius of curvature in polar coordinates. If 
the beam is centered and aligned with the optical axis we have 
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w(ijj) 
1 



7T 



V J (vo j» 

a 

WW 7 



|?/>(r)| 2 r 3 dr 
dip{r) 



Or 



dt/j*(r) 
dr 



r 2 dr 



7T 



/(VO ./o 



|^(p)|V^p 



(2.56) 
(2.57) 
(2.58) 



where = I(ip) is the integrated intensity in the transversal plane and p is the 
radial polar coordinate in the spatial frequency plane of the 2D Fourier transform of the 



field. It can be proved, using the integral transformation (2.37), that the three trans 



formation formulas are identical to the orthogonal case, and therefore the conservation 
(M 2 ) and the ABCD law remain valid for cylindrical symmetric beams. Once we have 
introduced this moment approach for the propagation of a general paraxial beam, it is 
useful to characterize the pointing stability, [28], of a laser beam with a misalignment 
factor |r/ m | 2 , which provides a global comparison between the misaligned beams and is 
expressed in the form 



!^^{x)[ip{x - S)e 



-ikax] * 



dx 



(2.59) 



I-oo \^{x)\ 2 dx 

for the one dimensional case, where r/ m is called the misalignment superposition 
integral, ip(x) is the field of the perfectly aligned beam and its misaligned copy is 
described by ij)(x — 5)e~ lkcix where 5 and a represent the transverse and angular shifts, 
respectively, of the misaligned beam in a given plane. The level of misalignment can be 
quantified by the number |r/ m | 2 , which takes a maximum unitary value for a perfectly 
aligned beam and goes to zero as the misalignment is increased. This number has all 
the desired features. It takes into account both transverse and angular shifts as well as 
the width and the divergence of the beam. In essence, this number quantifies by how 
much the misaligned beam differs in phase and intensity from the perfectly aligned 
beam in any plane, since it is invariant under propagation through an ideal optical 



system. For small shifts equation (2.59) reduces to 



1 - (M 



2\2 



e 2 + w 2 



(2.60) 
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where M 2 ,wq and 9q are the beam propagation factor, waist width and far-field 
divergence angle of the beam, in accordance with the second-order moments defini- 



tion (2.50), (2.42) and (2.46) respectively. In the derivation of (2.60) it was assumed 
that the transverse shift was measured at the waist plane. The variation range of the 
misalignment factor is < \r] m \ 2 < 1- The larger value of the parameters multiply- 
ing the transverse and angular shifts, means that the beam is more sensitive to the 



misalignment. It is important to stress that the formula (2.60) is a useful criterion to 



characterize the misalignment sensitivity of different beam geometry and this aspect 



will be investigated in Sec. 3.4.5 



2.8 Optical cavities with spherical mirrors 

An optical cavity Q or optical resonator is an arrangement of mirrors that forms a 
standing wave cavity resonator for light waves. Light confined in a resonator will 
reflect multiple times from the mirrors, and due to the effects of interference, only 
certain patterns and frequencies of radiation will be sustained by the resonator, with 
the others being suppressed by destructive interference. In general, radiation patterns 
which are reproduced on every round-trip of the light through the resonator are the 
most stable, and these are known as the modes of the resonator. Cavity modes are self- 
consistent field distributions of light, more precisely, electric field distributions which 
are self-reproducing (apart from a possible loss of power) in each round trip. The most 
common types of optical cavities consist of two facing plane (flat) or spherical mirrors. 



Now we have seen in Sec. 2.4 that the wavefronts of a Gaussian beam (and its higher 



order partners HG and LG beams) are paraboloidal surfaces with radii of curvature 



given by (2.27). Suppose that we fit a pair of curved mirrors to this beam at any two 
points along the beam in such a way that the radii of curvature of the mirrors are 
exactly matched to the wavefront radii of the gaussian beam at those two points. If the 
transverse size of the mirrors is substantially larger than the gaussian spot size of the 
beam, each of these mirrors will essentially reflect the gaussian beam exactly back on 
itself, with exactly reversed wavefront curvature and direction. These two mirrors thus 

1 Optical resonators are often called cavities. This term has been taken over from microwave 
technology (masers), where resonators really look like closed cavities, while optical resonators 
normally have a rather "open" kind of setup which doesn't really look like a cavity. 
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2.8 Optical cavities with spherical mirrors 



form an optical resonator which can support both the lowest-order gaussian mode, and 
(at slightly different wavelengths) the higher-order HG or LG modes as resonant modes 
of the cavity. Given two curved mirrors with radii of curvature R\ and R 2 , separated 
by a distance L, the calculation of the beam parameters follows directly from imposing 
that the wavefront curvature must match the mirror curvature at each mirror position. 
It is useful to introduce the resonator g parameters, g\ and g 2 , which are a standard 
notation in the field of optical resonators 

1 = 1-A and 02 = 1- A (2.61) 
rtl 112 

In terms of these parameters we can find that the self-reproducing gaussian beam 
has a waist size wq given by 



-o 2 = -\ , ^ ~° ig2 \ 2 , (2-62) 
and the spot sizes w\ and w 2 at the mirrors surfaces 

2 LX I g~ 2 2 LX I 9i / 0ftQ ^ 

w 1 = — ,/ — r and w 2 = — W — r (2.63) 

7T V ~ 9192) ~ 7T V g 2 (l - gm) 

It is obvious from these expressions that real and finite solutions for the gaussian beam 
spot sizes can exist only if the g±, g 2 parameters are confined to a stability range defined 
by 

< 9l g 2 < 1 (2.64) 

This condition is called stability range because this also exactly describes the condition 
required for two mirrors with radii R\ and R 2 and spacing L to form a stable periodic 
focusing system for optical rays. The stability of an optical cavity can be analyzed 
from the ray-optics point of view. Light rays that bounce back and forth between the 
spherical mirrors of a laser resonator experience a periodic focusing action. The effect 
on the rays is the same as in a periodic sequence of lenses or, for more general type 
of resonator, of ABCD optical elements. It can be easily proved that the rays passing 
through these equivalent periodic systems, are periodically refocused if the elements of 
the ABCD matrix obeys the inequality 

-1 < < 1 (2.65) 
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2.8 Optical cavities with spherical mirrors 



For a laser resonator with spherical mirrors the stability condition (2.65 ) reduces to 



( |2.64 ). The simplest resonator configuration to analyze are the symmetric resonators 



(baseline for Advanced LIGO) which have equal mirror curvature R and hence g pa- 
rameters <7i = g2 = g = 1 — L/R. The waist of the gaussian resonant mode is then 
obviously in the center of the resonator, with waist and end mirror spot sizes given by 



2 LX l+g 2 2 LX / 1 

Wq = h\w^) Wi = W2 = ^yT^7 ( 66) 

Another very often used configuration is the half-symmetric resonator (i.e. Virgo 
long arms cavities) in which one mirror is planar R\ = oo so that g\ = 1, and the other 
curved. The waist in this situation will be located on mirror number 1, with spot sizes 
given by 



" „..2 LA / 92 , o _ LX / 1 



Wn = Wi — J- and w 2 = — \ — y- r (2.67) 

° 1 vr V 1-52 vr y g 2 (l - g 2 ) 

The central point in the stability diagram, g\ = g 2 = 0, corresponds to the sym- 
metric confoca|^] resonator where the radii of curvature of the mirrors are exactly equal 
to the separation length L. The beam waist and spot sizes at the end mirrors are given 
by ^ ^ 

Wq = — and w\ = w\ = — (2.68) 

27T 7T 

The confocal resonator has overall the smallest average spot diameter along its length 
of any stable resonator, although other resonator may have smaller waist size at one 
plane within the resonator. Moreover the confocal resonator is also highly insensitive 
to misalignment of either mirror. 

We have seen that if a gaussian beam is found resonant for a given (infinite) spher- 
ical resonator, then all its Hermite-gaussian (or Laguerre-gaussian) functions can be 
resonant modes for that cavity. However for each of the transverse mode patterns, 
there are only certain optical frequencies for which the optical phase is self-consistently 
reproduced after each round trip (i.e., the round-trip phase shift is an integer multiple 
of 27r). These are called the mode frequencies or resonance frequencies and are equi dis- 
tantly spaced. The cavity modes are labelled by three indices in which the first refers 

1 This is referred to as confocal resonator because the focal points of the two mirrors coincide 
with each other at the center of the resonator. 
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to the axial mode number and the other two to the transverse mode indices. Using 
the expression for the total phase shift of a travelling gaussian beam, and imposing 
the boundary condition on the wavefronts at the mirror positions, we have that the 
resonance frequencies of the axial-plus-transverse modes in the cavity are given by 



where the + sign applies in the upper right quadrant (51,52 > 0) of the stability 
diagram and the — sign applies in the lower left quadrant. The frequency spacing, 
of the cavity axial modes is called free spectral range (FSR). The confocal resonator 
represent a situation where all the even-symmetry transverse modes of the cavity are 
exactly degenerate at the axial mode frequencies of the laser, and all the odd-symmetry 
modes are exactly degenerate at the half-FSR positions midway between the axial 
mode locations. This is the reason way, despite the optimum performance in term 
of misalignment stability, this configuration is not employed in GW interferometric 
detector, which must be able to operate on the single fundamental gaussian mode. 

2.9 Optical cavities: General Theory 

In the preceding section aperture diffraction effects due to the finite size of the mirrors 
were neglected. There, it was mentioned that resonators used in Fabry Perot optical 
cavities usually take the form of an open structure consisting of a pair of mirrors facing 
each other. Such a structure with finite mirror apertures is intrinsically lossy and, unless 
energy is supplied to it continuously, the electromagnetic field in it will decay. In this 
case a mode of the resonator is a slowly decaying field configuration whose relative 
distribution does not change with time [29]. The problem of the open resonator is a 
difficult one and a rigorous solution is yet to be found. However, if certain simplifying 
assumptions are made, the problem becomes tractable and physically meaningful results 
can be obtained. The simplifying assumptions involve essentially the quasi-optic nature 
of the problem; specifically, they are 1) that the dimensions of the resonator are large 
compared to the wavelength and 2) that the field in the resonator is substantially 
transverse electromagnetic (TEM). 




for HG 



for LG 



(2.69) 



(2.70) 
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So long as those assumptions are valid, the Fresnel-Kirchhoff formulation of Huygens 



principle can be invoked to transform the paraxial Helmholtz equation (2.4) with the 
appropriate boundary conditions for the cavity, into the familiar linear homogeneous 
Fredholm integral equation for defining the eigenmodes of an optical resonator. At a 
reference plane along the optical axis we have 



j n u n {x,y) = / K(x,y;x',y')u n (x',y')dx'dy' (2.71) 



where K(x, y; x' , y') is the Huygens kernel that propagates the complex optical wave 
amplitude u(x, y) through one complete round trip around the resonator. The kernel 
K(x, y; x' , y') in the Fresnel approximation can be given in terms of the ABCD matrix 



elements using equations (2.36) and (2.37) for cartesian and cylindrical coordinates 
respectively. 

The eigenfunctions u n (x, y) that satisfy this integral equation form the lowest-order 
and higher-order transverse modes of the resonator. Each eigenfunction u n for Eq. 



(2.71) has its corresponding complex eigenvalue 7 n . The fractional power loss per 
round trip for the n-th eigenmode is given by 1 — |7 n | 2 and the phase shift accumulated 
during a round-trip is given by arg(7 n ). 

Questions have been raised as to whether the lossy open resonator can be viewed 
as having resonant modes in either a physical or mathematical sense. The kernel of the 
Huygens integral above is symmetric but not unitary due to the boundary conditions 
over the finite mirror surface. Therefore, the different eigenfunctions u n are not power 
orthogonal in the usual sense. The orthogonality relations over the cross section are of 
the form f u n u m dxdy oc 5 nm rather than J u? n u m dxdy oc S nm . Hence, the total power 
J \ u\ 2 dxdy cannot be expressed as the sum of the powers carried by the individual 
modes. But, the eigenfunctions u n , although perhaps not normal modes in the usual 
sense, are still self-reproducing mode distributions in the sense that each such mode 
when excited will retain the same transverse form, attenuated only in over-all ampli- 
tude, after any number of transits around the resonator. Further, it seems clear both 
theoretically and experimentally that, as long as power is provided to the resonator to 
compensate for the losses, at least the lowest-loss mode, and in principle any of the 
higher-loss modes, can maintain a steady-state oscillation. In this sense it would seem 
that the functions u n are clearly meaningful transverse modes for an open optical res- 
onator. For optical resonators having rectangular mirrors (i.e., mirrors having different 
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radii of curvature in the x — y plane) the eigenvalue equation can be separated into two 
separate equations, each corresponding to the eigenvalue equation for a cylindrical strip 
mirror in the x or y variable. For mirrors of circular symmetry in the x — y plane the 
equation can be converted into cylindrical r, eft coordinates and then solved separately 
for modes of different azimuthal symmetry characterized by different exp(±im(f>) vari- 
ations. As a representative example we will recall the eigenvalue equation for optical 
cavities with spherical mirror which is of large use in GW community. In the following, 
and in future Chapters, we may refer to field amplitude at the mirror surface and not 
at a plane orthogonal to the optical axis^j For the case of circular mirrors, the equation 
(2.71) is reduced to the one-dimensional form by using cylindrical coordinate^] 



IpmRpmin) = - (j) \^ e - 2lkL £ riX^rur'JR^ir'Jdri (2.72) 



X m {n,ri) = I r 2 dr 2 J m ( — — )J m {— — ) exp 



(SlTi +gir 1 + 2g 2 r 2 ) 



where g\ and g 2 are the cavity's g-parameters defined in (2.61). 

This equation can be analytically solved only in the very special case gi = g 2 = 0, 
(a confocal resonator). In general must be solved by numerical analysis, as we will see 
in Sec. |3.2[ In the case of infinite mirrors the propagation between mirror surfaces is 
mathematically equivalent to a Fractional Fourier Transform, and the eigenfunctions 
of this linear operator are exactly the Hermite Gauss or Laguerre Gauss functions 



introduced in Sec. 2.6 The Fractional Fourier Transform has recently attracted much 



attention in optics studies as a tool to analyze paraxial optical systems imaging [30 . 

2.10 Optical response of a Fabry-Perot cavity 

For completeness we recall the basics principles of the Fabry-Perot optical cavities 
|31j(they constitute the fundamental configuration for the present and future GW in- 
terferometers and the focus of the R&D activity related to non Gaussian optics that 
-'in the paraxial approximation the transition from the field referred to the mirror surface 



to the field referred to a plane tangent to the mirror surface is given (as explained on pg 53 1 



by a phase transformation u{x,y) = v(x,y)e lk 2R t where R is the radius of curvature of the 
mirror and v and u are the field at the mirror surface and at the tangent plane respectively. 

2 The angular integrations can be computed analytically and led to the Bessel functions of 
the first kind and mth order. 
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will be analyzed in the following chapter). For simplicity we assume that the system 
is perfectly mode matched, that means the transverse field propagation law doesn't 
change inside and outside the cavity. Let consider the input field as a generic V>m 
and rj, tj and L are the field reflectivity, transmittivity and losses for each mirror [see 



Fig. 2.4 . The loss coefficient comprehends the absorption in the coating and any other 
source of optical power loss for the considered mode (r 2 + t 2 = 1 — £). 

¥trans 



¥ 2 



¥ 4 



¥3 



Figure 2.4: Electromagnetic fields in a Fabry-Perot resonator. 



The equations for the steady-state fields are: 



(2.73) 



l^trans 
^3 
I Aef 

where <j> is the phase gained (or loss) by the field in transit from mirror 1 to mirror 
2 (it includes the axial phase shift kL and the Guoy phase shift). The resulting field 



-^2^2; 

-ritpin + hlp4,. 



arc 



■01 

Iptrans 
tpref 



tie** 



1 — r^e* 2 ^ 

ht 2 e^ 
1 — ri^e* 2 ^ 



-n + 



t(r 2 e 



i(f> 



1 — ri^e 12 ^ 

The circulating power, P c i rc , stored in the cavity is given by 



(2.74) 
(2.75) 
(2.76) 
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t 2 

Pcirc =|V>i| 2 = t, r~27~T Pin (2-77) 



(1 - rir 2 ) 2 + 4nr 2 sin 2 (0) 
= 7 ^ n— = 9 2 PinA(<P) (2.78) 

where we have introduced the Fabry-Perot gain g and the Airy function A{4>). The 
maximum stored power corresponds to the peak of the Airy function, i.e. when the 
detuning phase 4> = nir, where n is an integer. The detuning phase can be changed 
to match the resonance condition both by varying the length of the cavity and/or the 
frequency of the input light. The half- width 5(f) of the resonances is determined by the 
finesse, F, of the cavity 

' with (2.79) 
2J 1 — r\r2 

The finesse is also related to the effective number of round-trips of the light inside the 
FP cavity by 3" « 7rN e f f and to the storage time of the cavity T s t Q by T sto 



■J2L 

'r*~i . 

7T C 
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Chapter 3 



Analytic and Numeric 
investigation of Optical Cavities 



3.1 Duality Relation 
3.1.1 Introduction 

We discuss here a duality relation between deviations from sphericity (or planarity) of 
mirrors in FP cavities. This duality proved a key point for our theoretical and exper- 
imental studies of cavities with non gaussian beam profiles [32 . In his recent work 
on a tilt instability for advanced LIGO interferometers [33 , P. Savov discovered nu- 
merically a unique duality relation between the eigenspectra of paraxial optical cavities 
with non-spherical mirrors: he found a one-to-one mapping between eigenstates and 
eigenvalues of cavities deviating from flat mirrors by h(r) and cavities deviating from 
concentric mirrors by —h(r), where h need not be a small perturbation. In the following 
section, we analytically prove and generalize this result. In this work, we prove this 
remarkable correspondence analytically, for an even broader category of cavities: those 
whose mirror shapes remain invariant under the parity operation, identified as spatial 
reflection in the two dimensional r-space (which is also equivalent to a 180° rotation 
around the cavity axis). Eigenmodes of such cavities can be put into eigenstates of 
parity, and we show that all corresponding eigenmodes of dual cavities have the same 
intensity profiles at the mirrors, with their eigenvalues satisfying 




(3.1) 
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where (— l) Pk is the parity of the kth eigenmode; subscripts c and f denote nearly 
concentric and nearly flat mirrors, respectively. We then illustrate its application to 
interferometric gravitational-wave detectors; in particular, we employ it to confirm the 
numerical results of Savov and Vyatchanin for the impact of optical-pressure torques 
on LIGO's Fabry- Perot arm cavities (i.e. the tilt instability), when the mirrors are 
designed to support beams with rather flat intensity profiles over the mirror surfaces. 
This unique mapping might be very useful in future studies of alternative optical designs 
for LIGO interferometers, when an important feature is the intensity distribution on 
the cavity optics. While such a duality relation is well-known between cavities with 
spherical mirrors, i.e., those with h(f) ~ ar 2 (for example see |24( [34"1 [35]). to our 
best knowledge no such relations had been established between generic cavities. In 
the following sections we present the work done in collaboration with E. D'Ambrosio. 
Another interesting approach to the duality relation was formulated by P. Savov and 
Y.Chen and is presented together with ours in a joint paper [32 . 

3.1.2 Analytical proof for mirror-to-mirror propagation 
3.1.2.1 Cartesian Coordinates 

In this section we focus on field distributions on mirror surfaces, and restrict ourselves 
to cavities with two identical mirrors facing each other. The extension to the non sym- 
metric cavity is presented later. We adopt the Fresnel-Kirchhoff diffraction formula to 
propagate fields from mirror surface to mirror surface (see e.g. [29] ). In this formalism, 
the field amplitude v± (r ') on the surface of mirror 1 propagates into 

v 2 (r) = fdV'X(r»i(r') (3.2) 



on mirror 2, via the propagator 

X(f, r ') = -^(1 + cos 6)e~ ikp k = ^, (3.3) 

Airp A 

from r' (on mirror 1) to r (on mirror 2), where p denotes the (3-D) spatial distance 
between these two points and 9 stands for the angle between the cavity axis and the 
reference straight line, as is illustrated in Fig. 3.1(a)| We know that the Fresnel- 
Kirchhoff integral eigenequation 

jv{r) = I dV X(r,f')v{r') (3.4) 
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h( F) 



h( F) 




(a) Nearly flat mirrors 




(b) Nearly concentric mirrors 



Figure 3.1: Geometric constructions for the propagators (3.6) and (3.8) 



univocally determines the eigenmodes v and eigenvalues 7 of the cavity. 
Applying the paraxial approximation 



0, p^L + 



2L 



h(f) — h(f') , 



and we can use 



ik 



-ikL ikh(r) 



e 2L 



\r-r'\ 2 ikh{r') 



(3.5) 



(3.6) 



in the integral eigenequation. 

Here the mirror surfaces deviate by h{r) from a flat reference, and the subscript f is 
used to reflect this convention. From here on, we will also refer to 3Cf as the nearly flat 



propagator. We now consider two slightly deformed concentric mirrors (see Fig. 3.1(b) ) 
so that the mirrors height with respect to the flat reference surface is 



h(r) = f 2 /L + b(f) , 



(3-7) 



where the height b{r) is the deviation from the concentric spherical surface (note that 
concentric spherical mirrors have their radii of curvature equal to L/2, and thus surface 
height r 2 /L). Inserting Eq. (3.7) into Eq. (3.5), we obtain the propagator for a nearly- 
concentric cavity, 

ik 



2irL 

e ikb(r) e +^\r+r'\ 2 e ikb(r') 



(3.8) 
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We use the term nearly concentric propagator for % b c {r, f '). Although we use the terms 
nearly-flat and nearly-concentric, h and b are not required to be small; in fact, they can 
represent any deviation from perfectly flat and concentric spherical mirrors. 

Now let us consider mirrors that are then invariant under parity, i.e., those in which 
we also have 

h(r) = h(-r) , 6(f) = 6(-f) . (3.9) 



so that 3Cf iC are both invariant under a spatial reflection 

(fjf*'} «-> {— r, — f '} 



and therefore, we have 



where we have defined 



(f) = v{—r) . 



(3.10) 

(3.11) 
(3.12) 



for two dimensional reflection. Equation (3.11) implies that all eigenmodes can be put 



into forms with definite parity. We derive the following relation between nearly flat 
and nearly concentric propagators, as constructed: 



that is equivalent to: 



nrh 

J\_£ 



-e 2ikL %- h {r,r') 



-e 2ikL X7 h . 



(3.13) 



(3.14) 



Suppose we have an eigenstate Vf of 3Cf , i.e., an eigenstate of a cavity with mirror 
deviating by (+/i) from flat surface, and we compute its eigenvalue j{ and know the 
parity eigenvalue (— l) p : 



9v f 



7f Vi , 
(-l)^f- 



(3.15) 
(3.16) 



(3.17) 



By applying Eqs. ( |3.14[ )-( |3.16[ ), we derive the correspondance 

which identifies v c = as the corresponding eigenstate of 3C~ ft , that is eigenstate of the 
corresponding resonator we denote the dual. The eigenvalue is 7 C ~- 



-2ikL 



• ,L (-l) p+1 7 f *. 
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Nearly Flat 


Nearly Concentric 




Kernel 


nrh 


<\ — h 


Ei| 


^enstate 


Vf 


V* 




Parity 


("1) P 


(-1)P 


Half-trip ei£ 


jenvalue 


7f 




Round-trip ei£ 


, r envalue 


m 


e~ AikL rf 



Table 3.1: Correspondence of propagation kernels, eigenstates, parities, and eigenvalues 
between dual configurations. 

We also induce that the parity is still (— l) p . The reverse is straightforward and the 
result is an established one-to-one correspondence between dual cavities. We summarize 



this mapping in Table 3.1 It is obvious to note that the corresponding eigenstates, Vf 
and , have the same intensity profiles on the mirror surfaces; for infinite mirrors, we 
know Vf(r) is real- valued (see Appendix b of [32]), so it is an eigenstate of the dual 
configuration itself. 

For cavities with identical mirrors facing each other, the full, round-trip propagator 



is just the square of the half-trip one. From Eqs. (3.13) and (3.11 ), we have 

2" 



rw-h 

J\_£ 



iikL 



-I, 



(3.18) 



which means that the same duality correspondence exists between eigenstates of the 
full propagator, with their eigenvalues related by 



Vc — e ?7 f 



(3.19) 



Note that when h(f) = r 2 /(2L) the two dual cavities are identical to each other. 
Using the relation that links the eigenvalues of two dual resonators, we can determine 
the spectrum 

7c = ± e -2^ 7 * = 7/ = e -ikL+inn/2 

where n £ N. The resulting separation between the eigenvalues is the Gouy phase 

e ie ° = e iaiccos { l - L /R) R = L 



computed for confocal resonators |24( 135] . 
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3.1.2.2 Cylindrically symmetric mirrors 



In most LIGO applications, cavity mirrors still have cylindrical shapes: h(f) = h(\f\). 
This allows us to decouple radial and azimuthal degrees of freedom, and simplify the 
eigenvalue problem. We shall follow roughly the notation of |35j. 
We adopt the cylindrical coordinate system: 



r(cos ip, sin if) . 



(3.20) 



Since % is now invariant under rotation along the z-axis, all eigenmodes can be put 
into eigenstates of rotation: 



v(r, (p) = R{r)e imLp , m = integer . 



(3.21) 



Inserting this into the eigenequation (3.42) and performing analytically the angular 



integration we obtain the radial eigenequation 



InmRnmir) = j K^^r, r')R nm (r')r' dr' , 



(3.22) 



where for each angular mode number m we have indexed the radial eigenstates by n, 
and 



K\ m) {r,r' 



i m+1 k j ( krr' \ j k \-L+h(r)+h(r')-^' 2 



(3.23) 



L V L 

is a symmetric radial propagator, in the nearly-flat description] Since K^ m) (r, r') is 
symmetric, we obtain orthogonality relations between radial eigenfunctions: 



f 

Jo 



(3.24) 



Using Eq. (3.7) again, for a configuration with b(r) correction from concentric spher- 
ical mirrors, we obtain the radial kernel of the nearly- concentric description: 



K b c(m) (r,r') = ^-J m l — \e 



ifc-L+6(r)+6(r') + ! 4f-' 



Comparing Eqs. (3.25) and (3.23), we obtain: 



(-1) 



m+l 



A f(m) 



e 2ikLj£-h 



c(m) 



(3.25) 



(3.26) 



1 f 27T 

1 Here we have used J n {z) = / e tzcostp e mLp dip, where J n {z) is the nth order Bessel 

2m 11 J 



function of the first kind. 



40 



3.1 Duality Relation 



This is a radial version of Eq. (3.14); here we know explicitly that all m-eigenstates 
have parity (— l) m . 

Following a similar reasoning as done in the previous section, for each angular mode 
number m, we can establish a one-to-one correspondence between radial eigenstates of 
a nearly-flat configuration to those of the dual configuration: 

[Rnm] c = [Rnm\*i ■ (3.27) 

The mapping of the eigenvalues are given by 

hnml = {-l) m+l e~ 2lkL [ 7nm ] f * . (3.28) 

Similarly, the round-trip eigenstates have the same correspondence, their eigenvalues 
related by 

[qn m ] c = e~ 4ikL [r? nm ]* . (3.29) 
3.1.2.3 Duality relation for non-identical mirrors 

In this section we will study the duality relation when the mirrors shapes are not iden- 
tical, but each still symmetric under a 180° rotation around the cavity axis. Since now 
the field distributions of eigenstates over the two mirror surfaces are not the same, 
we have to study the eigenvalue problem associated with the round-trip propagator, 
instead of the individual surface-to-surface ones. Nevertheless, we can still use the 



propagators (3.6) and (3.8) to build a system of integral equations relating field dis- 
tributions Ui(ri) and v 2 {r 2 ) over the two mirror surfaces. [Throughout this section, 
we use the subscripts 1 and 2 to refer to quantities associated with mirrors 1 and 2, 
respectively] If the mirrors deviate from parallel planes by h± 2(^1 we have: 



f A 2 r 2 % 12 {r ll r 2 )v 2 {r 2 ), (3.30) 
J Si 



7iui(ri 

>s 2 

72^2(^2)=/ dVi % 2 \{f 2 ,fi)vi{ri) , (3.31) 

JS! 

where 71 2 give the attenuation and phase shift experienced by the optical field in 
transit from one mirror to the other and 
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9Ci 2 (ri,r 2 ) 
X 21 (f 2 , n) 



ike 



-ikL 



2-kL 
ike~ lkL 
2-kL 



,ikh! (n ) - J| I r 1 - r 2 \ 2 +ikh 2 (r 2 ) 



,ikh 2 (r 2 )- # |r 2 -ri | 2 +ifc/ii (ri) 



(3.32) 
(3.33) 



are the propagators from mirror 2 to mirror 1, and from mirror 1 to mirror 2, 



respectively. The equations (3.30) and (3.31) give the field at each mirror in terms of 
the reflected field at the other but they can be combined to form the round-trip equation 
which states that the field at each mirror must reproduce itself after one round-trip. In 
the following, we will add a subscript f or c to make a distinction between quantities 
related to the nearly-flat or nearly-concentric case. 

h\h 2{ 



r) f vif(rx) = d r [ X^ f 2 (fi, r [) v lf (r [), 



(3.34) 

rifv 2f (r 2 ) = I d 2 r' 2 X h 2 p(r 2 ,r' 2 )v 2f (r' 2 ), (3.35) 
J s' 2 

where the common eigenvalue r/f is given by 71/72/ and the round-trip propagators 

/ d 2 f 2 %u f{n, r 2 ) X 2 i f (f 2 , ri) 
JS 2 



■Ap{r2,r' 2 



(1 



(3.36) 



In the nearly-concentric configuration, using kernels of the form (3.8) for the prop- 
agation from one mirror to the other and combining them as done for the nearly-flat 
configuration, we obtain the following nearly-concentric round-trip equation for the field 
distribution over the mirror 1 (similar formulas for the mirror 2 with the substitution 
1 <-> 2 throughout the argument). 



VcVic(n)= ^fi^^b^j^R) (3.37) 
Js[ 

^(nrt) = ~j s fr 2 e- 2 ^^) * • (3.38) 

. e 5f;|n+r2| 2 + ^|r2+r;| 2 +^ 

where b\^ 2 are the mirrors deviations from concentric surfaces. Using the assumed 
symmetry properties of the mirrors, the propagators for the nearly-flat and nearly- 
concentric cavity fulfills this relation 
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X^-^in,^) = e- AikL [X h x f\-r x ,-r' x 



1/ 

- 4 ^[3C^ 2 (ri,ri)]* (3.39) 



Equation (3.39), together with Eqs. (3.36) and (3.37), provides us with a more 



general duality relation, for cavities with non-identical mirrors: as long as the corre- 
sponding mirrors of two cavities A and B satisfy 

h«A(f) = j^-h aB (f), a = 1, 2 , (3.40) 
the eigenstates and eigenvalues of the two cavities will be related by: 

v aA = v* aB , r }A = e- mL rf B , a = 1,2. (3.41) 

3.2 Eigenvalues and FEM analysis 

In the preceding section we have analyzed an important relation which is satisfied by 
optical cavities in terms of eigenvalues and eigenfunctions of certain integral operator. 
In real application to GW Interferometers it is necessary to evaluate the eigenvalues and 
optical modes in situation where the mirror are not spherical and/or the finite size of 
the mirror become important. A few examples will be given below. We developed our 
own calculus packages in order to compute the optical modes and frequency spectrum 



of a FP cavity with non-spherical mirrors, which will be described in Sec. 3.6 This 
package found valuable applications also in the analysis of parametric instabilitj|^37 
of advanced detectors, due to its capability of obtaining, in a very efficient and precise 

1 Knowledge of the diffraction losses for higher order modes is important for predicting 
opto-acoustic parametric instabilities. These arise due to the resonant scattering of the cavity 
fundamental mode ujq against test mass acoustic modes u>t into other optical cavity modes ui\, 
which satisfy the condition lu ~ u> t + oji. The parametric gain R for this process determines 
whether the system is stable (Ro < 1) or unstable (i?o > 1)- The parametric gain depends 
linearly on the quality factor with higher order modes, which, on the optical side of the inter- 
action, depends inversely on the diffraction loss [36] - Thus it is very important to be able to 
obtain accurate estimates of modal diffraction losses in the optical cavities. 

2 In term of computational time. Other approaches like FFT propagation could obtain the 
same results but requiring up to two orders of magnitude more time than the approach described 
here (depending on the finesse of the cavity). 
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way, the eigenspectrum of an optical cavity. The equation which is the object of the 
next sections is 

jv(f) = / d 2 r'X(r, f')v(f') (3.42) 

J Mirror 

where the 2-D vectors r and r ' are the coordinates of the projections of mirror-surface 
points on planes orthogonal to the cavity axis at the mirrors positions, % is the round 
trip (or half round trip for symmetric cavity as Ad-LIGO) propagator and 7 and v(r) 
are the eigenvalues and eigenmodes (field distribution over the mirror surface) of the 
optical cavity. If the cavity mirrors are not equals, there is one eigenvalue equation for 
each field on the mirror surface; for the first mirror we have 



jvxin) = [ dViXf^ri^ri) (3.43) 
X^in,^) = - [ dV 2 e~ 2lkL ( -V) (3.44) 



5 2 



• exp |^-|fi + r 2 | - ^jj^2 + ^il + ihh\{ri) +ikhi(r 1 ) + 2%kh2{f-i, 

where h\ and /12 are the heights of the two mirrors respect with two planes orthog- 
onal to the optical axis and separated by a distance L. The equation for the transverse 
mode at the second surface is obtained with the substitution (1 *-* 2). These equations 
can be solved analytically only in very special cases (i.e. confocal spherical resonators, 
that is gi = gi = 0, [38]). In general some approximated solutions are found using 
different approaches. The finite-element method (FEM) has become a very powerful 
tool for the approximate solution of boundary-value problems governing diverse phys- 
ical phenomena. We developed some general numerical tools for the FEM analysis of 
optical cavities of interest for the GW community working on the optical set-ups. 

3.2.1 Cylindrical symmetric cavity 

We will now solve the problem stated above for a cylindrical symmetric cavity, which 
is of interest for the application on the ideal perfect optical cavities of the GW interfer- 
ometers both in their present design with spherical mirrors and in new configurations 
recently explored. It is obvious that for standard spherical optics the height function is 
given by h(r) = r 2 /(2R oc ), where R oc is the radius of curvature of the mirror. As long 
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as cavity mirrors posses cylindrical symmetry, we can separate the radial and azimuthal 
degrees of freedom, and simplify the search for eigenmodes. We adopt the cylindrical 
coordinate system: r = r(cos tp, sin <p). Since 3C is invariant under rotation along the 
cavity axis, all eigenmodes can be written in the form of 

v(r, (p) = R{r)e~ imLp , m = integer . 



Inserting this into the equation Eq. (3.43), and performing analytically the angular 



integration Q we have the reduced radial eigen-equation (a is the mirror radius) 

Tpmfli.pmM = / ' K^ ) {r 1 ,r J 1 )Ri,pm(r' i y i dr' 1 , (3.45) 



• exp | (rf + 2r\ + rf) + ikh x {rx) + ikhi(r[) + 2ikh 2 (r 2 ) i 

where for each angular number m we have indexed the eigenstates by p which is 
the radial mode number. 

The problem is thus reduced to a series of one-dimensional integral equation, one 
for each m. These are homogeneous Fredholm equations of the second kind and there 
are well documented standard techniques for their numerical solutions |39( 14*0] . 



There are two general methods of solving (3.45). One is based on iterative tech- 
niques applied to a discrete grid that extract the dominant or several of the lowest-loss 
modes; the classical Fox and Li algorithm falls in this category. 

For stable resonators the iterative method usually converges slowly because even 
higher-order mode losses are quite low. The second method, called Nystrom method, 
requires the choice of some approximate quadrature rules in order to reduce the integral 
equation to a matrix eigenvalue problem of dimension equal to the number of integration 
points N. The matrix method has two main advantages: it extracts the lowest N modes 
and their eigenvalues at one time, and its accuracy is determined by the order N of 
the matrix. It is certainly possible to use low-order quadrature rules like trapezoidal, 

1 f 2 * 

1 Here we have used J m {z) = — — / e lzcosip e" n(p dip, in which J m (z) is the mth order 

27ri" J 

Bessel function of the first kind. 
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mid-point or Simpson's rules, but since we are looking for a quite accurate solution of 
the eigenvalue problem (the modulus of the eigenvalues must be calculated with high 
precision for a realistic estimation of the diffraction losses) this would require a large 
number of integration points and therefore a lot of computation time. We used the 
Gaussian quadrature rule in order to have a greater accuracy with fewer points. 

Typically, the value of the radial coordinate is bound by the radius of the mirror. 
We denote as a this maximum value. Let x% be the abscissas and Wi the weight factors 
for GaussLegendre quadrature in the range (0,a). This technique has fine meshing at 
the edges of the computational window and more course meshing in the center. 



The equation (3.45) for a certain m can be approximated as 



N 

'yR(xi) = K(xi, Xj)xjWjR(xj) (3-47) 
Defining Ri the vector R{xi) and K{j the matrix K(x{, Xj)xjWj we have transformed 



the integral eigenvalue problem Eq. (3.45) to a matrix eigenvalue problem r yR = KR 
which can be solved for 7 and R by any standard mathematical package. The con- 
vergence of the eigenvalues and eigenvectors of the discretized problem to those of the 



continuous eigenvalue equation (3.45) for sufficient large N is guaranteed by mathe- 
matical theorems |41j . 

Once the eigenvectors are determined, a particular mode may be evaluated at any 
point r (not a sample point) by calculating the sunj^J 

N 

r,Xj)xjWjR pm (xj) (3.48) 

3 

The complex eigenvectors vi- pm (r,(p) = Ri- P m{f)e~ tmip correspond to the resonator 
mode patterns at the first mirror surface. The fractional power loss per transit of the 
mode due to diffraction effects at the mirrors is given by £ = 1 — |7 pm | 2 . As is usual in 
optical resonators, the eigenvalues can be sorted in decreasing magnitude order. The 
phase shift (p pm is given by the angle of 7 pm , which is the phase shift suffered (or 
enjoyed) by the wave field in a roundtrip. The resonant condition requires that the 
total phase shift along the axis of the cavity be an entire multiple of 2ir rad; thus, 

1 This formula has been demonstrated to be the best interpolating function between the grid 
points. This is called Nystrom interpolation formula. 
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separating the longitudinal phase shift^J (f> pm = —2kL + (3 pm = —2irq. The resonance 
frequencies of the axial-plus-transverse modes in the cavity is given by 

Um = k{q+^) (3.49) 

where /o = ^ is the free spectral range of the cavity, i.e the frequency spacing be- 
tween successive longitudinal resonances. We implemented the described procedure in 
a Mathematica® program, using its built-in routines for the solution of eigenvalue 
problerrj^J 

In Fig. |3.2| we show that the CPU time for the computation of ten eigenvalues as 
function of the grid size N. The advantage of using the Gaussian quadrature rule is 
evident if we consider that a 1000 grid nodes in an equal spaced mid-point rectangular 
quadrature rule can achieve an accuracy of 10~ 6 in the magnitude of the eigenvalues; 
by comparison a Gaussian quadrature rule with 100 nodes reaches an accuracy of 10~ 13 
in the eigenvalue magnitude value. 

3.2.2 Advantages with respect to iterative method 

The advantages of the eigenvector method are: 

• The propagation matrix has to be calculated only once, then the main task is to 
solve the set of eigenvalues of the transit matrix, and each eigenvalue can directly 
derive one set of eigenvectors, which just represents one mode distribution on the 
mirror, without hundreds of iteration. 

• Phase shift and amplitude loss per roundtrip can be easily evaluated from the 
eigenvalues. 

• The multi-modes can be obtained with a single calculation. 

• There is no need of an initial field distribution, which has to be carefully given 
in the FoxLi or in the Prony method to get the convergent results. 

1 The minus signs are due to the sign convention in the propagator. 

2 As a technical remark is useful to point to the fact that when diffraction losses are very 
small, i.e. when the absolute value of the eigenvalues is close to one, it may be necessary to in- 
crease the significant digits in all the computations and one of the advantage in Mathematica® 
is that it can handle approximate real numbers with any number of digits. 
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Figure 3.2: CPU time for the computation of 10 eigenvalues for Ad-LIGO like cavity. 
[PC Laptop ,CPU 2 GHz, RAM 780 MB] 

As an example of the capabilities of this simulation tool let us mention the use 
in the calculation of realistic diffraction losses in the Ad-LIGO FP cavities, for the 
estimation of parametric instability. In eleven lines of Mathematica® code we are able 
to compute the diffraction losses of ten higher order modes in less than four seconds of 
CPU time using 100 grid nodes. Barriga et al. |37j used an iterative process based on 
FFT algorithm, which requires one houiQof CPU time for the computation of each mode 
diffraction loss. As it is shown in Fig. |3.3[ the agreement between the two calculations 
is very good with a maximum deviation of around 5% due to intrinsic limitation of the 
FFT algorithm for this kind of analysis. 

3.2.3 Arbitrary mirror shape 



In the preceding paragraph we analyzed the problem (3.43) in the very special case 
of cylindrically symmetric mirrors and ideal mirrors with no imperfection. In real 
application, we have to analyze broader situations, where the mirrors deviate from 
Barriga in private communication. 
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Diffraction Losses Comparison for LG Modes 



1 8 
Optical modes 



■ Eigen Values * FFT a Clipping 



Figure 3.3: Comparison between diffraction losses calculated using eigenvalue FEM 
algorithm of this thesis, FFT program and clipping approximation in the case of Ad- 
LIGO FP cavities with mirrors diameter of 31.4 cm. [Courtesy of P. Barriga] 

the ideal shape by some amount, i.e. intrinsic manufacturing defects or deformations 



due to thermal effects. It is thus very useful to have an algorithm to solve Eq. (3.43) 
without the restriction of Sec. 13.2.11 In literature there are two well known technics 
for approximate solution of this 2D problem both based on iterative procedure; to find 
the steady-state mode distribution, an arbitrary initial wave distribution is propagated 
repeatedly around the resonator, undergoing changes from transit to transit and loosing 
energy by diffraction. After many iterations, the wave converges to a steady-state 
stationary mode. In one case the propagation between mirrors was carried out by 
expanding the optical wave in an eigenfunction expansion using as a basis set the 
Hermite-Gaussian beam functions characteristic of free-space optical beam propagation. 
Because these functions are normal modes for free-space propagation, this basis set has 
both conceptual and computational advantages, as demonstrated in [42 . However, this 
expansion has also the same practical disadvantage characteristic of most such normal 
mode expansions. If one expands a discrete N x N point two-dimensional function in 
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an M x M series of basis functions, the amount of computation required goes up in 
proportion to (NM) 2 , or in proportion to iV 4 if M = N. The other approach, |43| . is 
based on a Fast Fourier Transform algorithm to compute the propagation between the 
surfaces and will be described in more details in the next Section. The FFT method has 
the advantage of requiring a smaller computational effort which increase as N 2 log 2 N 
respect with the iV 4 of the other iterative approach in the real space. Nevertheless if 
the diffraction losses are small, the iterative process could converge very slowly and 
require a long computational time. For this reason we developed a procedure for the 



calculation of the eigenmodes, that reduce the problem (3.43) to a matrix eigenvalue 



problem which can be solved by standard numerical routines. The problem is that in 



a simple discretized version of (3.43) over two 2D grids, one for each mirror, we have 



an equation between 2D matrix related by a 4D kernel. We use the following trick to 
reduce it to a standard matrix problem. First we choose a numerical integration rule 
for the evaluation of integrals in a compact bi-dimensional domain (mirror surface). 
We will use, for simplicity, a mid-point rectangle quadrature formula^ such that the 
mirror is sampled at equally spaced points in the x, y coordinates and assume that the 
mirror have equal radii o. We first select the dimension of the 2D grid N 2 and the 
spacing between nodes in this grid 



A x = A. 



2a 



(3.50) 



y ~ N 

It is however a waste of memory to keep track of the node points that fall outside 
the mirror boundary; we therefore build a M x 2, with M < N x N matrix, that 
contains only the coordinates of cells inside the mirrors 



\XM VM/ 



(3.51) 



The equation (3.43) is then reduced to the matrix equation 

M 

iv k = ^2 K kj v v k = l---M 



(3.52) 



1 The extension to more efficient quadrature formulas in polygonal domain is quite straight- 
forward. 
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where the matrix K is given by 



/ , n 2 M 

Ku =-{ — - I A 2 e - 2ikL \^ e -^hV. x k'n) 2 +{yk-yi) 2 +^i-^j) 2 +(yi-yj) 2 

k3 \2ttLJ x ^ 
v 7 i=i 

ik [hi (x k ,y k )+h 1 (xj ,yj)+2h 2 (x t ,y t )] 



(3.53) 



This eigenvalue equation can be solved by any standard mathematical package. We 
implemented the above described routine in a Matlab® code because of its efficiency 
in dealing with large matrices. With respect to the ID grid of dimension N, the size 
of the matrix K increases as iV 4 , requiring a notable amount of memory to store its 
values and compute the eigenvectors. 

3.2.4 Constraints on FEM analysis 



It is interesting to point out some constraint of FEM analysis applied to Eq.(3.43) 
since it is often an overlooked aspect in papers related to this technique for optical 
simulation. The grid size N is set by a compromise between accuracy of the calculations 
and computational effort required. The kernel in the integral equation is a rapidly 
oscillating function of the coordinates and sampling it at the grid nodes introduces 
inevitably some errors. A faithful discrete representation of this function is obtained 
applying a sort of sampling theorem of signal analysis. The fundamental constraint is 
that the phase difference between two adjacent point must be less than ir. There are 
two contribution to this phase, one independent on the mirror profiles and the other 
directly given by the height of the mirrors. They give the following constraints on the 
maximum spacing allowed between grid nodes 



LX 



and Ah [Ax] < - (3.54) 

max 2 



These constraints set, for a given quadrature rule, the number of minimum grid 
points; it is important to point out that the first constraint is determined by the laser 
wavelength and the cavity length, whereas the second constraint is set by the mirror 
geometry and the laser wavelength. Once these minimum requirements are fulfilled, a 
further increase in iV results in a better accuracy in the extraction of the eigenvalues 
and eigenvectors. 
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3.3 FFT 

A new package based on Mathematica® was developed to simulate an optical cavity 
when realistic mirrors are included. It can moreover deal with many features includ- 
ing mirror tilts and shifts, beam mismatch and/or misalignment, diffractive loss from 
finite mirror apertures, mirror surface figure and substrate inhomogeneity profiles, in 
particular, using deformation phase maps that are adapted from measurements of real 
mirror surfaces and substrates, as well as fluctuations of reflection and transmission 
intensity across the mirror profiles. Similar programs [44; 45 developed by many other 
research groups, have been used for a variety of applications by the gravitational wave 
community, including numerous design help and performance estimation tasks. We 
decided to build our own code because we needed a fast and reliable tool for the op- 
tical simulation of our mesa beam cavity; the time we would have spent in learning 
and modifying the existent programs was employed to write and debug our own code. 
Incidentally, thanks to its simplicity and versatility, the programs revealed itself very 
useful for investigations related to Ad-LIGO and Virgo thermal compensation system. 

The required procedure is to expand the propagating wave at any given plane into a 
series of plane-waves or spatial- frequency components, travelling at different angles in 
k— space. The component plane waves are propagated from one transverse plane to the 



next by a simple propagation formula, Eq. (2.12). The transformations from coordinate 



space to k— space and back again are formally equivalent to Fourier transformations, 



Eq. (2.10 ), and can be very efficiently carried out using the fast Fourier transform (FFT) 
algorithm. Fourier transformation carried out using the fast Fourier transform (FFT) 
algorithm provides the opportunity for a significant increase in computational efficiency 
The computational labor required to expand a two-dimensional N x N function in a 
two-dimensional N x N Fourier expansion only increases as N 2 log2N, which is a much 
smaller number than N 4 for moderate-to-large ./V values. However, propagation from 
plane to plane is not sufficient for our purpose. In general we need to represent the 
propagation from the surface of one mirror to the other, which could be of standard 



spherical shape, or even non-spherical as we will see in Sec. 3.6. One fundamental point 
for the application of this algorithm is that the action of a mirror on the optical field, 
can be modelled without using a diffraction integral. Following jH] it is convenient to 
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discuss this issue in the Fourier space. If ip(x, y, 0) is the field in the plane z = 0, and 
ijj(x, y, z) the field on the mirror surface, we have 

i>{k x , k y , z) = e~ iz V k2 - k i- k2 y tp(k x , k y , 0) (3.55) 

The argument of the imaginary exponent can be expanded keeping in mind that we 
are dealing with paraxial optics 

/ z(kl + k 2 ) 

zy/k*-k*-k**kz- K x 2k v > (3.56) 

In order to estimate the order of magnitude of the second term with respect to the 
first one let us remind that the divergence of a beam is defined as the second moment 



of the Fourier transformed field, Eq.(3.68); Therefore we have 



< k l + k D , el zM 2 

wvr — ~ = (3-57) 

2k Z7T ZrTT 

where we have used the definition of the generalized Rayleigh range zr introduced 



in Sec. 2.7 Since the second term is a factor 9 2 smaller than the first one, we can 



conclude that for surface height z <C zr, with good accuracy we can write 

${k x , k y , z) = e~ ikz 4>{k x , k y , 0) (3.58) 

The expression of the field on the mirror surface of equation z = f(x, y) is related 
to the field at a reference plane z = by 

Vw/(x, y) = i> P lane{x, y)e~ lkf ^ (3.59) 

with this approximation is obvious that the reflection and transmission operators 
from a mirror with non planar surfaces are described by 

R{x, y) = re i9{x ' y) d{x, y) and T(x, y) = te ie ^ v) d{x, y) (3.60) 

where t and r are the ordinary scalar amplitude transmission and reflection coef- 
ficients and g{x, y) and e(x, y) are the functions representing the local phase change 
due to reflection or transmission trough a mirror (these functions depend on the mirror 
surface shape and on the material refraction index). 
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Therefore we can represent propagation in free space by a phase factor in the Fourier 



space according to equation (2.15), and the interaction with a mirror by phase factors 
in the direct space. 



A simple FP cavity like Fig. 2.4 can be thus modelled using just multiplication of 
elements of matrices and FFT routines. For example the intra-cavity field tpi obeys the 
steady-state equation 

^ l = T 1 i) in + R liT PR 2 ,iPe iLp i)i where P = $- 1 G L 9 and G L = e~ lkL e iL ^^- 

(3.61) 

where 7 indicate the Fourier transform, P is the propagator down the cavity length 
L, the reflection operators distinguish between right (r) or left (/) incidence and ip is a 
uniform phase added to achieve the fine tuning of the cavity (it plays the same role as 
a microscopic length change). Let us assume that we choose a square grid in the real 
space of dimension N x N and the corresponding grid in the Fourier space. Thus each 
mirror reflection or transmission operation is reduced to a pixel-by-pixel multiplication 
between a mirror operator map and the e- field slice on the reference plane near to it; the 
propagation between two mirrors is given by an FFT, followed by a element-by-element 
multiplication between the transformed matrix and the Gl matrix, and then an inverse 
FFT. The simulation program models a static optical cavity, which means that we do 
not follow the time evolution of the field inside the cavity but we iterate the equation 



(3.61 ) until a steady-state is obtained for the field. The simple field relaxation algorithm 
requires a pre-specified threshold of accuracy: the iteration process ends when the field 
calculated at the M th iteration is equal, within the accuracy requested, to the field 
calculated at the M — 1 th iteration. Other more sophisticated relaxation algorithm are 
described in [45 j and will not be adopted in this work. 



3.3.1 FFT limitations 

The propagation calculation using FFT procedure require proper attention to questions 
of sampling and aliasing. These problems have been deeply discussed in (35] and we 
will point the fundamental issues. Let a being the mirror radius and consider a square 
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grid of dimension N 2 and side length W which contains the mirror at its centeiQ Let 
Ax being the spacing, in the x/y directions, between the grid nodes (W = N Ax). 
These parameters set the k-space grid length k max and spacing Ak 



vr _ttN 2tt _ 2tt . _ N N 

Ax W ' N Ax W ' 2 ' ' 2 

(3.62) 

It is well known that sampling a function which has power at spatial frequency 
above the Nyquist frequency, = 1/(2 Ax), causes aliasing. One possible solution is 
increasing the highest spatial frequency in order to push this limit above the spectral 
distribution of the function. To quantify the requirements for the /c-space gridding, 
let us remember that k x /k gives the propagation direction of one Fourier component 
(plane wave) and that the beam divergence, Go is the second moment of the Fourier 
transformed field. In paraxial beam propagation, the k-spectrum of the beam is there- 
fore limited by 7O0, where 7 is a constant of order unity. The maximum propagation 
angle and the ^-resolution pose the following constraints to the grid choice 

kmax ^ ... NX 780 A 
— — >7@o^W< — — , (3.63 
k 2760 ^ 7@o 

Remembering that for a generic paraxial beam, M 2 = WqQqtt / \, we see that the 
second constraint is easily satisfied, but the first needs some care: it requires a large 
number of grid nodes N 3> 1 but also constrains the resolution Ax of the grid by 
Ax <C It is moreover important to consider the aliasing effect due to large- angle 

scatter from mirror deformations: a very simple way to take care of the power scattered 
by one fourier cell into the next replica, is to increase the window size W well above 
the mirror diameter, and apply the so called "zero-padding" method. It simply forces 
to zero the power outside the mirror surface. We thus have that for most application, 
W = 4a is a good choice. There is also another aspect that has to be considered among 
this program constraints. Since the actions of the mirror are represented by pixel-per- 
pixel phase multiplication, the optical path difference between two neighboring pixels 
must be smaller than A/2. This sets some limitations on the mirror shape and tilt 
we are able to simulate with a given grid finesse. For our optical cavity (L = 7.32m, 

1 For the moment we do not set a constraint on the relative size W/a. It will be a crucial 
point in the aliasing discussion later. 
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N = 128, Ax = 0.35 mm, a = 1.3cm) we have the limitation 1.5 mrad on tilt we can 
model. 

3.4 Mesa beam 
3.4.1 Introduction 

The fundamental Gaussian TEMqq has great flexibility, is used in many applications, 
and has an easy and attractive mathematical form. Nevertheless, many cases require 
different intensity profiles: a cylindrical volume with a diameter d = 3w is just a quarter 
filled by a gaussian intensity distribution, therefore a gaussian beam is not appropriate 
for high power laser applications. The case of interest here is that of the reduction of 
thermal noise of mirrors. Most gravitational wave interferometric detectors (GWID) 
measure the variation in phase between light beams resonating in two perpendicular 
cavities caused by the passage of a gravitational wave. Any physical displacement 
of the reflective surfaces of the cavity mirrors also creates phase variations and thus 
contributes noise to the measurement. In particular, random displacements of the test 
masses' reflective surfaces due to ther mo dynamical fluctuations is a major source of 
fundamental noise in the frequency range of maximum sensitivity for next generation 



interferometers (Fig. 1.3). 



Gaussian laser beams are typically used to measure the position of the mirrors in 
these new detectors. Because of the peaked nature of the gaussian beams, they are 
not well suited to average over the localized thermal fluctuations. It is possible to 
significantly reduce the test mass thermal noise using modified optics ( "graded-phase 
mirrors [IS]") that reshape the beam from the conventional gaussian intensity profile 
into a flat-top beam profile. Similarly, flat topped beams would be of interest to reduce 
mirror thermal noise in metrology and in frequency and distance standards. 

One set of potential beam profiles that satisfy the request of a flat-top, more uniform 
transverse power distributions, is given by the supergaussian functions, with analytical 
form SG(r) = e~ rl . Such a field pattern can be obtained shaping a gaussian beam with 
a diffractive optic, and use a refractive lens to produce a focused flat top intensity. In the 
case of gravitational wave interferometers, the optimized displacement sensitivity also 
requires very small diffraction losses (~ 10 ppm per bounce) in their Fabry-Perot arms. 
This precludes ordinary refractive techniques to flatten the resonant beam shape in the 
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cavities. It was shown that certain optical resonators formed by two facing graded phase 
mirrors - i.e. aspherical mirrors - can have a fundamental mode with supergaussian 
features [47J. The mirror profile, in this case, can be designed to match the phase 
wavefront for a prescribed intensity profile, while maintaining small diffraction losses. 
A different analytical approach to design flattened beam was proposed by Gori [48j: 
the flat top beam can be obtained as a sum of Laguerre Gauss modes in the cylindrical 
coordinate system. The main advantage of such transverse field distribution is that the 
field that they produce upon propagation can be evaluated in a simple way without 
introducing any approximation except, of course, that the paraxial regime is assumed 
to hold. Tovar [49] has proposed a new class of beams, called multi- Gaussian beams, 
that can reproduce the flatness of a super-gaussian transverse field. They consists of 
a small sum of finite-width gaussian beams side by side with the same width, phase 
curvature and absolute phase. Unlike the flattened Gaussian beams, each of the multi- 
Gaussian beam components can be traced individually without resorting to further 
series expansion. Hence, this approach has an analytical form more desirable to study 
the beam diffraction characteristics. 

Furthermore, beam propagation must be exactly known in order to characterize 
interferometer performances. 

Thorne and others proposed and theoretically studied [E2 [20] [50] the possibility 
of using a particular class of aspherical mirrors, the so-called "Mexican-hat" (MH) 
mirrors, in a Fabry-Perot cavity to generate a wider, flat top laser beam, the mesa 
beam. Such a beam is predicted to significantly reduce all sources of test mass thermal 
noise by better averaging over its surface fluctuations pJ2 |20j [50] EU [52] [53]. A 



complete analysis of this problem will be given in Sec. 4.8 



3.4.2 Mesa Beam and Mexican Hat mirrors cavity 

The mesa-beam must have an intensity distribution that is nearly flat across the light 
beam, and then falls as rapidly as possible (to minimize diffraction losses) at the beam's 
edges. For a given length L, and laser wavelength A there is a particular Gaussian beam, 
called minimal- Gaussian beam, whose radius increases by a factor \pl between the beam 
waist (at the cavity's center plane) and the cavity's mirrors, which provides a minimally 
spreading Gaussian beam. The mesa electromagnetic field, in a quasi-flat cavity, is a 
superposition of minimal- Gaussian fields whose axes are parallel to the cavity axis and 
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lie within a cylinder of radius b centered on the cavity axis. The (non-normalized) field 
distribution at the mirror plane is 



U FM (r) 



dVe~ 



<b 



2ir I c 
'o 



, (r*+r")(l-i) 
2 ™0 



rr'(l 



r' dr' (3.64) 



where Io is the modified Bessel function of zeroth order and wq = s/L/k is the 
waist of the minimal-Gaussian (where k is the wavenumber). For a cavity to have 
the mesa beam as an eigenmode, the surfaces of the mirrors must coincide the mesa 
field's surfaces of constant phase at that distance from the waist (this is true for infinite 
mirrors, but also a quite good approximation for finite mirrors if the diffraction losses 
are small) . The resulting height distribution as a function of radius r is given by 



hMH(r) 



Arg[U{r)\ - Arg[U(0)] 



(3.65) 



The resulting mirror's radial profile is shown in Fig. 3.4(b) along with the spherical 
(nearly-flat) mirror profile which support a Gaussian beam with the same diffraction 
losses; notice the shallow bump in the middle and the flaring outer edges which resemble 
a Mexican hat (sombrero) and gave the mirror its name. 



Power 




height 



Mirror profile 



r 

w 



1 2 3 4 5 6 

(a) Intensity distributions (a.u.) for a Gaus- 
sian beam (red) and a Flat Top beam 
(blue) having the same diff. losses. 




(b) Mirror profiles: spherical mirror (red) 
supporting the Gaussian beam and Mexi- 
can Hat mirror (blue)supporting the mesa 
be am. 



Figure 3.4: Comparison between a Gaussian beam and a mesa beam. 

The corresponding optical design has shown strong tilt instability at high power [33 than 
the corresponding Gaussian beam. K.S.Thorne proposed a different version of the mesa 
beam, that is supported by nearly concentric and opportunely shaped mirrors; this new 
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version provides the same intensity profile at the cavity mirrors surface (and thus the 
same thermal noise), but introduces a weaker tilt instability (better than Gaussian 
beam cavities with the corresponding nearly concentric spherical mirrors analyzed by 
Sigg and Sidles [54 ) — as calculated by Savov and Vyatchanin [33] and by us in 
Sec. |3.5.2 A general method to design a family of optical cavities has been proposed 
by Bondarescu and Thorne [SS], from nearly flat resonators to nearly concentric ones. 
This kind of mesa beam is constructed by coherently overlapping minimal Gaussian 
beams with axes in diverging directions sharing a common point at the cavity waist. 
The non-normalized field distribution at the waist plane (cavity middle point) is given 
by 



j r waist 
U CM 



1 



r dr 



2tt 



+ 2irpr cos(0) 

>e w o 



w 



r -2 rb 



(3.66) 
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This field can then be propagated along the optical axis using the Huygens prop- 



agator (2.37) adapted to a simple free space propagation. Examples of normalized 



power distributions of nearly flat and nearly concentric mesa beams are plotted in the 



upper panels of Fig. |3.5| In these plots, we take b = 4wo, which corresponds to the 
configuration proposed for Advanced LIGO . 

Even if the beam profile doesn't conserve the same shape along propagation (see 



Fig, 3.4.2 and 3.4.2), which is the case of Gaussian beam, we will show in the next 



section that is possible to characterize the beam trough the moments definitions of 
Sec.l27n 



3.4.3 Analytical investigations of mesa beams 



Mesa beam have been studied mostly numerically [20J and Sec. 3.6.1.2 perturba- 

It 



tively [13 ED] and Sec. or using an expansion on Gauss Laguerre functions 
is however of great utility to have some analytical formula for studying the beam 
characteristics and beam propagation trough an optical system without recurring to 



numerical analysis. In this section we follow the guidelines given in Sec. 2.7 for the 



definition of some parameters which are the usually measurable quantities (2.56) in an 



optical beam and which behave in a simple way under a generic ABCD transformation. 



The expression of the mesa field as given by equation (3.64) seems quite intractable 
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from an analytical point of view, in particular the integral cannot explicitly computed. 
However it is a convolution product of a gaussian with the disc function of radius b and 
in the Fourier domain this become a simple product of the transformed functions, which 
are well known. The first important step is to express the formula for the generalized 
beam width and radius of curvature in the Fourier domain. Since we are dealing with 
cylindrically symmetric problems it is useful to work with these definitions of Fourier 
transform and Hankel transform 



f(k x , k y ) = dxdy e'^+^fix, y) = f(s) = 2tt / rdr f(r)J (rs) = H(f) 
Jr 2 Jo 



(3.67) 



y^ = 4^J R2 dk * dk v e l{kxX+kyV) f(k x , k y ) = f(r) = i 



2tt 



sds f(s)J (rs) = H~\f) 



We recall here the formula of the generalized beam parameter assuming a normalized 
field t/> [] 



W = 2\ TT 
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dr 
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(3.68) 
(3.69) 
(3.70) 



Using the definition (3.67) and the properties of Bessel functions is easy to prove 
the following equalities 



2vr 



2vr 



r 2 ^(r)Jo(rs)rdr = 

r — - — jQ\rs)rdr 
or 



d 2 ip(s) 1 dtp(s 



+ 



s ds 



ds 2 ' s ds 



It is then possible to express the equations (3.68) in a different way 



(3.71) 



1 Note that the coefficient of the divergence contains different it factors due to the different 
definition of the Fourier transform. 
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(3.72) 
(3.73) 
(3.74) 



First of all we need to find a normalization for our flat mesa (FM) beam (3.64) 



otherwise the above formulas are not analytically computable. We use the Parseval 
relation^ applied to the transformation (3.67) in polar coordinates 



poo 1 poo 

dxdy\f(x,y)\ 2 = 2TT r dr\f(r)\ 2 = — sds\f(s)f 
R 2 JO 27T JO 



(3.75) 



Lets consider the mesa beam (3.64) at the beam waist 
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(3.76) 



Now, the Hankel transformation of the gaussian and of the disc function are well 
known 



H(e ™o) 
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and H(D(b)) = 2nb 



Ji(bs) 



(3.77) 



The integral of the intensity in the Fourier domain can be performed analytically 
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(3.78) 



We have therefore the normalized field at the beam waist Vo fm 



WfibVYTTZ JD(b) 



^o,FM(r) = - ^_ 3 I e ™o dx dy Q 



b 



_ , °e H k- x2) iJ xdx (3.79) 

lr The same technique has been used in |57j . Our results confirms their calculation of the nor- 
malization factor and correct some errors for the Gaussian beam coupling efficiency calculated 
there. 
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The expression of the Fourier transformed field is 

^o.fm(s) 



2vr ™fc 2 J] (6a) 
e 4 — ^ — L 



(3.80) 



The expression of the FM beam in other plane along z is just a straightforward 
application of the propagation rule for each Gaussian composing the beam. The calcu- 
lation led the following resultj^] 



ipz,FM.(r, z) 
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where zq is the Rayleigh range and w z is the beam radius at the z plane of the 
minimal Gaussian field 
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(3.82) 



The Fourier transformed field is 



~ 20™ o / .z \ Ji{bs) 
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(3.83) 



Using the expression (3.80) in equation (3.72) we can calculate the FM beam waist 



Wq,fm an d the divergence 0q,fm- 
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(3.84) 



(3.85) 



Using (3.83) with (3.72) we can calculate, at any plane z, the beam width and the 



radius of curvature 



1 Omitting the phase factor common to all the minimal Gaussian given by the Gouy phase 



shift 



— z arctan — 
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It is possible to verify that these quite complicated propagation rules can be put in 
the standard form 
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(3.88) 
(3.89) 



3.4.4 Concentric Mesa beam 



In the Concentric Mesa (CM) beam case, (3.66 ), we have the great advantage of having 



an explicit form of the field at the waist plane. The normalization is thus immediately 



performed using the same integral (3.78). 
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(3.90) 



and the propagation is obtained using the Huygens integral (2.37) with L = B 
z, A = D = 1, C = 0. 
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It interesting to observe that the field ipo,CM( r ) has exactly the same for of the 
Fourier transformed ^o,fm( s )- This is a manifestation of the duality relation between 
these two kind of beam; one beam is related to the other by a transformation which 
has the form of an Hankel transformation Q As seen for the FM case with equations 

1 Further details about this aspect of the duality relation are given in [32] 
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(3.88) and (3.89) it is sufficient to calculate the beam waist and divergence to obtain 



the propagation formulas. For the CM it is more useful to work with the field (3.90) 



and use the first of (3.68) and the last of (3.72) for this computation. The results are 
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(3.92) 
(3.93) 



3.4.5 M 2 parameters and misalignment sensitivity 



The parameter M 2 , introduced in (2.50), is an invariant of the beam along the prop- 



agation and is a very useful tool in beam propagation analysis. At the beam waist is 
given simply b}|]] 
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Combining (3.84) and (3.85) we have 
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In Fig, 3.9 is shown Mp M as a function of the integration disc radius b. 



Looking at the formulas for the CM beam (3.92 ) and (3.92 ) we obtain the important 
result 
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which is another manifestation of the duality relation. 

The misalignment sensitivity of a generic paraxial beam is quantified by the mis- 



alignment factor \rj m \ (a, 5) introduced in [28] and briefly recalled on pg. 26 (where 6 



and a represent the transverse and angular shifts, respectively, of the misaligned beam 
1 For a pure Gaussian beam M 2 = 1. 
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in a given plane). The misalignment superposition integral of Eq. 2.59 can not be an- 
alytically performed for FM and CM beams but the approximated expression Eq. 2.60 
offers a good approximation for small misalignments. 

We will therefore continue the misalignment analysis using the approximated ex- 
pression 

W ._ 1 _ (J0 »(^ + ^) (3.97) 

where Go and Wq are the beam divergence and beam waist respectively. This factor 
is trivially computable for Gaussian beam. To analyze the mesa beam misalignment 



sensitivity in both configuration, FM and CM, we plot in Fig. 3.11 the beam divergence 
and beam waist for different value of the constituent parameter b. 

The misalignment sensitivity ratios between the FM and CM configurations are 



represented in Fig. 3.12 



3.4.6 Gaussian beam coupling to FM and CM beams 

For the application to advanced GW interferometers FP cavities it is of fundamental 
importance to find the best Gaussian injection beam which has the maximum coupling 
with the Mesa beam resonating in the Mexican hat cavities. This can be computed 
analytically using the above formulas. Lets take the scalar product, at the waist plane, 



between a gaussian field of waist w„ and the FM beam field (3.76) 
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After the integration we have 



Qfm = ^Lw n II - e 



(3.99) 



Fig. 3.13 shows the power directly coupled, Q 2 FM from a gaussian beam into a FM 
beam. The maximum coupling, which occurs for w g solution of 
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(3.100) 
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in agreement with [20] , can be quite large, i.e for b = 4wo we have a power coupling 
of 94% with a gaussian beam of w g ~ 3.62u;o- 

We now proceed to the same calculation for the CM beam, which seems the most 
promising design for advanced interferometers 
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Carrying out the integration we have 
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As shown in Fig. 3. 14 for the proposed design 6 = 4u>o, a power coupling of 94% 
occurs with w g ~ 0.28wq. 

Therefore, direct coupling of a gaussian beam with a mesa beam can be very efficient 
in both configurations. 
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Figure 3.5: Comparison between nearly flat (left panels) and nearly concentric (right 
panels) Mesa beams. Upper panels: normalized intensity profiles at the center of the 
cavity. Middle panels: normalized intensity profiles at mirror surfaces Lower panels: 
phase fronts at the position of the mirrors. 
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Figure 3.6: FM beam profiles along propagation on free space, z is in unit of zq 
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Figure 3.7: CM beam profiles along propagation in free space. z is in unit of zq. The 
intensity distribution are scaled to the on-axis value 
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(a) Transverse shift. 





1 








° 


8 

6 


I 

I 

i 
i 


\ 1 ??m 1 2 

\ — " l^ml^pp 

A 


CM 




i 


l\ 


~i 
fr 


4 


/ / 
/ 1 
/ I 


1 \ 
1 \ 





2 


/ i 
/ i 
/ i 


1 \ 
I \ 
1 \ 
1 \ 



L_i , ^L. _1_ , I , _L ^5* L. 

-1 -0.5 0.5 1 

a/e 

(b) Angular shift. 



Figure 3.10: Comparison between the exact and approximated expression of \rj m \ 2 (a, 5) 
for FM beam with b = 4wq . 
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(a) Divergence. 
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(b) Waist width. 

Figure 3.11: The scaling to the divergence and waist of the minimal Gaussian reveals 
the duality relation between the FM and CM design. 
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Figure 3.12: Ratios between translation and rotation sensitivities for FM and CM. 




Figure 3.13: Power coupling between a Gaussian beam and a FM beam. 
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Figure 3.14: Power coupling between a Gaussian beam and a CM beam. 
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3.5 Application of the duality relation to real problems 

Spherical cavities are not optimal in terms of their thermal noise: (the two types of) 
mesa beams, whose intensity profiles are natter given the same loss specification, turn 
out to provide much lower thermal noises. For these beams, the larger the parameter b, 
the lower the thermal noises, but the higher the diffraction loss. The loss specification 



of Advanced LIGO corresponds to b = 4wq which is the case we study in Fig. 3.4 
While having the same diffraction losses and thermal noises, dual configurations do 
differ significantly in a very important aspect: their eigenspectra are different. Thus, 
any problem using modal analysis of optical cavities will reveal these differences. 

It has been recently pointed out [54] that the force induced by radiation pressure, 
when the end mirrors of an optical resonator are misaligned in tilt, can induce a positive 
feedback and further increase the tilt. This effect is obviously power dependent and 
increases with power. In the case of anti-symmetric tilt, the coupling between the 
laser radiation pressure and the misalignment makes the mirrors tilt further, while in 
the symmetric case the torque induced by the light beam impinging on the mirrors, 
counteracts the tilt. Calculations show that the antisymmetric tilt is much more critical, 
and is therefore advantageous of reducing the effects of the antisymmetric case even at 
the expense of worsening the symmetric case. Although control systems are designed to 
keep the mirrors aligned, they have limited authority, and minimized coupling between 
radiation pressure and cavity misalignment is desirable, when there is positive feedback. 
Quantitative assessments of the problem [54] (that depends on the geometry of the 
resonator) have shown that optical cavities with almost flat mirrors are more prone 
to the antisymmetric instability with respect to almost concentric cavities (where the 
comparison is done between geometrical configurations that correspond to the same 
transverse distribution of the fundamental mode on the reflective surface of the mirrors) . 



In Fig. 3.15 there is a geometrical explanation of this phenomenon for the spherical 
mirror case. 

Because of this issue, the current baseline design for Advanced LIGO has been 
changed from nearly flat to nearly concentric [58 . 

It is important to evaluate the radiation pressure problem for the case of mesa 
beams. We will use our proved duality- relation to show that, in the quasi-concentric 
Mexican-Hat configuration, the cavity is less prone to become unstable and therefore 
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easier to control [59]. For general, non-spherical cavities, a perturbative approach 
must be used to calculate the tilt instability. The torque due to radiation pressure is 
expressed in terms of numerically found eigenvalues and intensity profiles of the cav- 
ities 's spatial eigenmodes. Using the numerical work done for the mesa beam in the 
nearly flat configuration, and using the mapping of eigenvalues and eigenmodes between 
the nearly flat and nearly concentric cavities (to obtain the estimation of the tilt insta- 
bility for nearly concentric Mexican-Hat cavities without having to solve the eigenvalue 
problem again) we proved that, for the antisymmetric case, the nearly concentric MH 
configuration is much more stable than the nearly flat MH one, and even slightly more 
stable than the corresponding nearly concentric spherical configuration proposed for 
Advanced LIGO. The same results were obtained by Savov and Viatchanin [33] using 
a similar modal analysis. 

3.5.1 Perturbation theory for finite size mirrors 

Since the propagation between finite mirrors is described by non-unitary operator, 



Eq.(3.43) we can't apply the standard perturbation theory, based on a set of power- 
orthogonal modes. We derived a rigorous procedure for perturbative analysis in the 
finite mirror case, which converges to the standard approach in the limit of infinite 
mirrors. With this method we can apply a perturbative approach even in those cases 
where the diffraction losses are large and the optical modes are greatly affected by the 
mirror finiteness. The drawback is that, since the modes of an open resonator are the 
eigenmodes of a non-unitary operator, it cannot be rigorously guaranteed that they 
form a complete set. We simply make the assumption that these modes can be used 
as a basis set and check this assumption numerically by examining the convergence of 
the series expansion. 

In general we have that for a non hermitian operators, a biorthogonality relation 
holds between the eigenfunctions of the operator K = K(r, rb), the transposed operator 
K T = K(fo,r) and the hermitian adjoint operator K H = K*(fo,r). K,K T and K H 
will have separate and different sets of eigensolutions fulfilling these relations 



K<f> n = J n (f> n K T if n = K n (f n K U n = a n UJ n 

a* = K n = 7 n and ip n = io* n 
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but no simple relationship between eigenfunctions cp n and iv n for non-hermitian 
operator (K ^ K H ) 

Instead of being power-orthogonal (J <j)* n (r)(j) m {f)dr ^ S nm ) the eigenfunctions of K 
will be biorthogonal to the eigenfunctions of the corresponding transpose (or hermitian 
adjoint) operator 



(p„{r)fm(r)df= J (t> n (r)u)* m (r)dr = 5 r , 



In the special case of simple optical cavities the kernel (3.43) is always symmetric 
in r r 1 and this implies 4> n = tp^\ 

In the following we will use a compact formalism for hermitian scalar product and 
usual integration and we will assume that the eigenmodes are power-normalized even 
if not power-orthogonal. 

(u\u) = J dru*(r)u(r) = 1 (3.104) 
[u\v] = J dfu(r)v(r) (3.105) 
[u\G\v] = I ' drdr'u{r')G{r,?)v(r) (3.106) 



To keep light notation we will indicate simply with n the set of integers which 
characterize the eigenvalues and eigenvectors of our problem (e.g. p, m for unperturbed 
aligned cylindrically symmetric cavity). 

Assume that we solved the eigenvalues problem for the kernel Ko, which corresponds 
to the unperturbed cavity 

Ko^n = li°Un (3.107) 

and we want to study the effect of a small perturbation (e.g. small mirror tilting) 
so that the resulting full kernel can be written as 

K = Ko + K\ and the eigenvalues equation as (Ko + Ki)ip n = 5 n ip n (3.108) 

1 This is not true for more general types of resonator like ring resonator. These cases are 
discusses in I2H. 
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We'd like to solve this problem using a perturbation series and to keep track of 
powers of the perturbation we will make the substitution K\ — > eK\ where e is assumed 
to be a small parameter in which we are making the series expansion of our eigenvalues 
and eigenvectors. 



5n = ^ + e^ + e^ + .... (3.109) 
i> n = N(e) («/>„ + c nk(e)<f>k) (3.110) 

Cnk(e) = ecU + e 2 c% + .... (3.111) 

where the superscript (0),(1),(2) are the zeroth, first, and second order terms in 
the series and N is a normalization factor. 
Using these equations we can write 



{Ko + K^^ + Y.Cnkie)^) = (4° ) +^k 1) + ^ { n ) + ----)(^n + J2 C nk^)4>k) (3-112) 

kj^n kj^n 

For this equation to hold as we vary e, it must hold for each power of e. Looking at the 
first three terms (e , e l ,e ) 



K (/>n = l^Un (3.113) 

K E + ^ £ 4V, = 7i 0) E c i 2 >* + ^ E + i^k 

k^n kj^n kj^n k^n 

(3.114) 

The zero order term is the solution of the unperturbed problem and so there is 
no new information there, but we can extract the first and second order correction to 
eigenvalues and eigenvectors from the other two terms using the orthonormality of the 
functions fa. Multiplying these equations by <p n or <fo. and integrating over the free 
coordinates we have after straightforward calculations 
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7 




(i) 



[(f) n \Ki\^ n \ 

[<t>n\4>n] 
[<t>k\Ki\<f>n] 1 



7, 



4°)_ 7 (°) [<f>M 

i ^ a^ii<M) 2 



ttM&n 7^ "if 



(3.115) 



Where the denominators take care of the normalization properties for the eigenfunc- 
tions. The normalization factor N(e) played no role in the solutions to the eigenvalues 
equation since that equation is independent of normalization. We can easily derive the 
normalization factor for the first order corrected eigenfunction . 



The correction is of order e and can't be neglected at this level of approximation. 

These results are nearly formally the same as in quantum mechanics but we have to 
keep in mind the different meaning of our notations (the expression of the normalization 
constant is completely different due to the fact that we don't have power-orthogonal 
functions). 

3.5.2 Evaluating tilt instability 

The case of antisymmetric mirror tilt, as indicated above, is the possible source of 
laser-power driven instability. 

The mechanical torque exerted by the misaligned fundamental mode on the mirror 
is given by 




(3.116) 




(3.117) 
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For a symmetric spherical mirror cavity, the torque can be computed, at the leading 
order, using very simple ray-beam physics. The pure geometric contribution is given 



by T = 2P/c8x where 5x is the displacement of the optical axis (see Fig. 3.15 given by 
L0/(1 — g). This very quick result is completely in agreement with |33( 154] . 

In the special case of mirror's antisymmetrical tilt by a small angle 9 we have that 
the full kernel describing the propagation from one mirror's surface to the other is 



j^tilt (-f f ) = — c ik6r cos M+» fcfc ( r )- 4| \r-r'\ 2 +ikh(r')+ik0r' cos(y') (3 118) 

flat flat i 2irL 



If ^ ^ ~ta = 2™ we can wr ite this kernel as 

Kft* t -flat(r, r) = K (r, f>) (l + ikOr' cos(^) + ikOr cos(p)) + O{0 2 ) (3.119) 
and we can identify 

Kx(f, r) = K (f, r)ikd(r' cos(ip') + r cos(^)) (3.120) 

We now calculate the lowest order correction to the fundamental eigenmode and 
eigenvalue. 

7^ oc [</> |ifi|<fo] = [4>o(r)\K (f, r )ik6(r' cos(ip') + r cos(^)) \(f>o(r')] 

= [M^)\K (ry)ik9r' cos(^)l0o(^)] + [<fa(f)\Ka(f, 7)ikercos(<p)\<j> ^)] 

= 2-ff\k0 [(j)o{f)\r cos((p)\4>o(r)} = since 4>o(r) = (j) (r) 

so the first order correction to the lowest eigenvalue is zero and before calculating the 
second order correction let us calculate the first order coupling coefficient 



(1) _ [0fcjjfjJjV] 1 _ 

° ok ~ 7 f_ 7 (°) [4>k\4>k] ~ 

i 0) -l { k 0) [Mrl\Mrl] 
= 7r fl 7 fc 0) +7 o° } [Mr)\rcos(<p)\(l>o(r)] 

7 (o)_ 7 (o) [umm 
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Remembering the form of the unperturbed eigenfunctions <^ pm (r) = R pm {r) cos(nnp) 
we have that the possible non zero terms involve only <^> p i(the radial index may vary) 
and the result is 

(l) _ .^ Tp? +7qq [^ P i(r)|rcos(y?)|0oo(r)] /q 19 -n 

(2) 

Since the same kind of integral enters in the definition of 7o we can easily see that 
the sum over k ^ reduces to just one contribution 

(2) , 2fl2 (7 P ? +t£ } ) 2 ([0 pl (r-)|rcos(y)^oo(r)]) 2 
700 " 7 g)- 7 g) [0oo(r1|0oo(r1] 

We now calculate the first order corrections in the p = case (lowest radial eigen- 
functions) so that the eigenfunctions entering in the integrations are of the form 



4>oo{r) = R 00 {r) <j>oi(r) = Roi(r) cos(^) 
From Eq. |3.117 and the above results we have the dominant contribution is given 



by 



T ~ ^ 23R (c^ 1 (0oo|rcos(^)|0oi)) (3.123) 
This can be written as explicit function of the eigenvalues and radial eigenfunctions 



c V 700 -Toi J drr R 01 (r)R 01 (r) J ) c 

(3.124) 

From our numerical calculation for the solution of the eigenvalues equation using 



the algorithm explained in Sec. 3.2 follows these results: Fig. 3.16 shows the intensity 



profiles of the lowest order modes of a symmetric Mexican hat symmetric cavity (both 



FM and CM). Tab. 3.2 shows an example of eigenvalue calculation for a Mexican hat 
mirrors cavity like Ad-LIGO proposal (a = 15.7cm, b = 4 • wq) 

Using our equivalence relation we can immediately write the eigenvalues and eigen- 
vectors in the nearly concentric case 



CMH _ (_-i\m+l( FMHy -2ikL 
1pm \ I \ipm ) " 
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7pm FMH 


m = 


m = 1 


p = 


0.759684 - 0.650267? 


0.835831 - 0.548805i 


p = 1 


0.979863 - 0.199425?' 


0.984052 +0.175851? 



Table 3.2: First eigenvalues for a nearly fiat Mexican Hat cavity with L = 4Km, 
a = 15.7cm, A = 1.064/im and b = Awq. 



R™ H (r)=[R™ H (r)Y (3.125) 
and the numerical calculations for the overlapping integrals are performed using the 



eigenfunctions of the FMH configuration (Fig 3.16). 

We calculate the ratio between the coupling coefficients of the nearly- flat and nearly- 
concentric configurations 

O-CMH _^ 1 
OiFMH 247 

which means that the MH-concentric configuration is about 250 times less sensitive 
to antisymmetrical-tilt induced torque than the MH-flat configuration. 

In Fig. |3.17 we investigate the sensitivity of different mesa beam design with respect 



to the baseline Advanced LIGO configuration. 



As we will see in Sec. 4.8 increasing the mesa beam radius has the advantage 
of reducing the mirror thermal noise. In particular the configuration with b ~ 4.65wo 
should be the optimal choice for Ad-LIGO in terms of thermal noise, and this correspond 
to a factor of 470 between the torques in the FM and CM configuration. 
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Figure 3.15: Each of these optical resonators is formed by two identical spherical mir- 
rors, with radius of curvature much larger (nearly fiat cavity) and smaller (nearly 
concentric) then the distance between the optics. When the mirrors are tilted in a 
common mode (on the left) by an angle 6 the optical axis shifts by R6. The result- 
ing torque is smaller in nearly concentric cavities than in nearly fiat resonators. The 
torque restores the cavity alignment when the mirrors are tilted in a differential mode 
(on the right). The change in the optical axis is a rotation 6/(1 — ^) around the center 
of the cavity. The resonating beam is displaced by ^f/(l — 2r) on the optics. This 
displacement occurs in opposite directions and tends to restore the correct alignment. 
Since this rotation is larger in nearly concentric cavities, these benefit more from the 
restoring torque. Quantitative assessments can be found in literature and substantiate 
these considerations. 
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TEMqo TEM10 




Figure 3.16: First Mexican Hat cavity modes (power distributions in a.u.): they replace 
the Gauss-Laguerre modes of the spherical cavities. 
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Figure 3.17: Radiation pressure torque sensitivity ratios of flat Mexican hat (FMH) to 
concentric Mexican hat (CMH) mirrors cavity. 
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3.6 Mexican Hat cavity experiment 

This work is a collaboration with M. Tarallo , J. Miller, R. DeSalvo, E. D'Ambrosio, 
P. Willems, B. Simoni, J.M. Mackowski and A. Remillieux, colleagues at University of 
Pisa, the LIGO Laboratory at Caltech and LMA (Laboratoire des Materiaux Avances) 
in Lyon. The main motivation of this project |60( I6T] is to demonstrate the feasibility 
of building a Fabry-Perot cavity with Mexican Hat mirrors, which can support a Mesa 
beam. 

The aim of this experiment is to explore the main properties of a single optical cavity 
before any eventual use in a second generation gravitational-wave interferometer. In 
particular, we are interested in studying the experimental mesa field achievable with 
realistically imperfect mirrors and how its behavior differs from that of a gaussian 
field with respect to perturbations such as cavity misalignments. Other groups have 
demonstrated mesa beam cavities using deformable mirrors |62| [6"5] , but such mirrors 
are not obviously usable in low-noise gravitational wave interferometry. 

The MH mirror production technique sets the main constraint to the prototype 
geometry. 

The production of the Mexican hat mirror has been undertaken in the Laboratoire 
des Materiaux Avances in Lyone (LMA). They use a three steps deposition process 
over a micro-polished flat substrate: the general shape coating, the corrective coating 
and the multi-layer coating. In the first step a "rough" Mexican hat shape is deposited 
with a precision of about 60 nm using a profiled mask and rotating the substrate to 
generate the cylindrical symmetry. The mask, calculated from the thickness profile of 
the ideal Mexican hat, is placed between the sputtered flow of silica and the rotating 
substrate. 

The second step is a more precise correction of the "general shape" previously 
obtained. This method controls the deposited profile with a precision of about 10 nm 
Peak-to- Valley (PV). Nevertheless, it is not possible to coat more than 100 nm with 
this technique, because it would require a deposition time which is too long. The 
measurement of the achieved mirror shape after the first deposition is performed using 
an interferometric technique. The comparison between the achieved and the desired 
mirror shape generates a data file which is used to move the robot arm that positions 
the mirror in front of the corrective silica beam. The main limitations of the corrective 
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technique come from the measurement of the wavefront, the precision of the robot arm 
movement and the size and resolution of the SiC>2 corrective beam. The maximum 
achievable slope is 500 nm/mm, thus setting a limit on the smallest feasible Mexican 
hat mirror. Finally, a high reflectivity Si02/Ta20s multi-layer coating is deposited on 



the corrected substrate. The slope of the Mexican Hat mirror profile 3.4(b) sets the 
radius of the smallest feasible mesa beam made using our technique to about 6 mm. 

Our smallest practical mirror size sets our cavity length to ~ 16 meters. We further 
reduced the physical length of the structure to ~ 8 m by building a half-symmetric 
cavity (a single MH mirror paired with a flat mirror at what would be the midpoint of 
a full length cavity) and then to ~ 4 m by folding. In this way it was possible to build 
a rigid suspended cavity. 



Fig. 3.18(a) shows the suspended cavity. Three Invar rods fix the cavity length to 
L P rototype = 2x3.657 m = 7.32 m, with a folding mirror at one end of the structure and 
the input and end mirrors on the other end. Five triangular spacers maintain structural 
rigidity, with the outer two spacers bolted at the ends of the structure and containing 
the mirror mounts. The cavity is suspended by two pairs of mar aging steel wires 
from GAS (Geometric-Anti-Spring) [64 blades, providing both horizontal and vertical 
isolation. The whole is suspended in an aluminum chamber for thermal stability and 
protection from air currents; this chamber is not evacuated. 

The test MH mirror was designed using the waist size of the minimal gaussian with 
L = 2L pro totype as a reference length, so that the resulting mesa beam had a radius 
r mesa = 6.30 mm. We required that our test mirror have similar diffraction loss around 
its aperture as an Advanced LIGO test mass; in order to have 1 ppm diffraction loss 
the mirror radius was set to a = 13 mm. 

Due to the technical difficulties of the MH figure deposition on the flat substrate, 



the MH mirror had non-negligible figure error [see Fig. 3.6 , mostly in the central bump 
where the height is just 27 nm. In particular, the figure error reaches a maximum of 5 
nm at the edge of the central bump. 

Further details on the experimental set-up and the development of the project are 
given in the MS thesis of Simoni [65J and Tarallo [66 . 
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MH mirror 

(b) Schematic diagram of the cavity prototype 

Figure 3.18: Mesa beam cavity prototype. 
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Figure 3.19: Deviation from idea MH shape in the central mirror area, x-y pixel size 
is 0.35 mm. 



3.6.1 Comparison between experiment and simulations 

In order to simulate our FP cavity we used a dedicated version of the FFT implemen- 



tation described in Sec. 3.3 Its aim is to investigate the impact of misalignments and 
mode matching with the input Gaussian beam driving the FP cavity and to study the 
cavity behavior with imperfect mirrors using real mirror maps for setting requirements 
and tolerance for Mexican hat mirror manufacture and control constraints. From the 
mirror's maps we can make theoretical prediction of the sensitivity of this new type of 
optical cavity to mirror imperfections and alignment. Since the geometry of the Mex- 
ican hat mirror of this FP cavity is completely different from the standard spherical 
mirrors supporting Gaussian beams, the effects of misalignment for mesa beams are 
more important. In the Mexican hat case, the change in the surface of the mirror due 
to orientation has a more general effect: not only does the cavity have a new optical 
axis but the phase profile sensed by the beam is quite different. For spherical mirrors 
this is not the case since an incident beam experiences the same curvature at all points 
on the mirror surface. 
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We have mounted the sample 5008, one of the three made in the LMA laboratory. 



As shown in Fig. 3.6 it presents significant deviations from the design mirror shape. 
In particular, there is an evident slope in the central part of the mirror, which man- 
ifests itself by modifying the expected beam shape, as shown in the FFT simulation 



of Fig. 3.20(a) . We used the map of this sample mirror in our simulation program to 



characterize the properties of the resonating light. Our simulations have shown that 
the resulting beam shape could fit the expected flat top behavior: by slightly changing 
(~ rad) the alignment of the mirror we can partially correct for this intrinsic tilt. The 
best beam we could obtain with this imperfect mirror is shown in Fig. 3.20(b) [ assuming 
perfect flat input and folding mirrors). The ripples in the central area are inevitable 
due to the roughness of the mirror surface (limited by the accuracy of the corrective 
coating deposition) and set a limitation of about 10% PV on the flatness of the power 
distribution on the top of the beam. 



y(cm) 





(a) Before intrinsic tilt correction. (b) After mirror tilt of 1 /irad. 

Figure 3.20: FFT simulations with the real mirror map of the sample 5008. 
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3.6.1.1 Alignment and intrinsic defect correction 

The greatest experimental difficulty was found in obtaining a sufficiently precise align- 
ment to achieve a flat top power distribution in the cavity. In a cavity made with 
spherical mirrors, a tilted mirror presents the same spherical profile to the opposite 
mirror, but shifted sideways, and the cavity simply resonates the same transverse mode 
spectrum centered upon a shifted optical axis. In contrast, the MH mirror has a non- 
spherical shape, and any misalignment destroys the cylindrical symmetry of the cavity. 
In such a situation the resonant beam senses a mirror with a suboptimal profile and the 
cavity mode will thus have a radically different intensity distribution and phase front. 
When our cavity was in such a misaligned state, higher order and distorted modes were 
found easily. We were able to reproduce these patterns with our simulation program 
by changing the alignment of the MH mirror and/or shifting the input Gaussian beam 
as shown in Sec. 6.3 of 166 . 



Fig. 3.21 shows the simulated cavity spectrum when a tilt of 1 fii&d is applied to the 
MH mirror. The higher order modes are less effected by the MH mirror tilt with respect 
to the fundamental mode; this is traduced in the experiment in a greater difficulty in 
aligning the mirrors to obtain a fundamental flat mode. On the other side higher order 
(misaligned) modes are obtained very easily. 

3.6.1.2 Cavity modes and spectral distribution 

The extreme sensitivity to misalignment proved to be caused by our nominally flat 
folding and input mirrors, which had surface deviations of order 60 -j- 100 nm from 
flatness before we fixed aluminium compensating/strengthening rings to them. Having 
reduced the flat mirrors' flexure, it was possible to lock the cavity to a stable funda- 
mental mode with nearly uniform distribution. In fact, taking into account the residual 
warping of the flat mirrors and MH mirror imperfections, this mesa beam is consistent 



with the best achievable using our current prototype MH mirror [see Fig. 3.22 . Figure 
|3.23| shows four beam profiles. Two are experimental data, smoothed with a 0.1 mm 
(5 pixels) gaussian kernel to clean them of digitization noise and dust diffraction rings. 
The other two are profiles simulated using our FFT code. The simulated profiles repre- 
sent the leakage field at the output bench (~ 5 meters from the input mirror) achieved 
applying the ideal corrective tilt at the MH mirror. 
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Figure 3.21: Effect of 1 /^rad tilt on the simulated cavity mode spectrum. 



Analyzing the profile in Fig. 3.23 the normalized absolute power in the simulated 
profile not fitted by the mesa TEMoo is 3.4%. By comparison, the normalized absolute 
power in the experimental profile not fitted by the mesa TEMoo is 3.8%, after the 
normalization of the power on the CCD, with a peak to valley deviation from the flat 
profile of about 9.4%. These numbers suggest that the resulting mesa beam is very 
close to the experimental limit due to the imperfect MH test mirror. 

In Fig. 3.24| we show some higher-order transverse mesa beam modes observed. 



These modes are superficially quite similar to the Laguerre-gaussian modes for a spher- 
ical mirror Fabry-Perot cavity. However, in the precise power distribution there are 



differences as shown in Fig. 3.25. The experimental data is in good agreement with the 
expected mesa TEMio profile. As for the fundamental mode, there is some asymmetry 
due to the mirror imperfections. 
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100C- 




Figure 3.22: Three-dimensional profile of the mesa fundamental mode. Experimental 
data. 



Experimental beam profile. x-axIs We ' -jcperimental beam prufile : y- c 




Figure 3.23: One-dimensional profiles of fits to the mesa beam profiles. The top row 
shows normalized experimental data as measured at the CCD camera. The dashed line 
is the best fit mesa profile. The bottom row shows profiles extracted from the FFT 
simulation with best corrective tilt applied. In this case, the transverse scale is taken 
at the MH mirror. 



93 



3.6 Mexican Hat cavity experiment 




Figure 3.24: High order mesa beam transverse modes: (a) TEMio, (b) TEMn, (c) 
TEM20 • Experimental data. 
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Intensity (digits) 




— 1 — 1 — ' = x (mm) 

12 3 4 



Figure 3.25: The mesa TEMio profile (thick black line). The light gray line show 
the theoretical mesa TEMio, which better fits the data than a Laguerre-Gauss TEMio 
mode (dashed line). 

A simple study has been done to compare the numerical predictions for our res- 



onator and cavity length sweeps. Results are shown in Tab. 3.3 The second column 
is calculated assuming perfect MH and flat Input mirrors. Although the estimate has 
large uncertainty (0.2 MHz due to frequency spacing and systematic effects) we can 
conclude that the numerical predictions are well respected. 

The mirrors imperfections and the quite large deformations of the Input mirror cause 
a deformation in the mode shapes and some shifting in the resonances spectrum. The 



simulation shown in Fig. 3.26 is performed using our algorithm described in Sec. 3.2.3 
The cylindrical symmetry is broken by the mirror's imperfection/deformations and the 
power distribution over the lobes of the higher order modes is no more symmetric (as 
experimentally observed) 
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Peak 


Af exp (MRz) 


TEMp/ (expected) 


Af(TEM pl ) 


1 





00 





2 


0.4413 ±0.2 


01 


0.4141 


3 


1.198 ±0.2 


02 


1.0945 


4 


1.574 ±0.2 


10 


1.6542 


5 


2.144 ±0.2 


03 


1.9905 


6 


2.900 ±0.2 


11 


2.8789 


7 


4.161 ±0.2 


12 


4.1754 


8 


4.414 ±0.2 


20 


4.4050 


9 


5.549 ±0.2 


13 


5.5523 


10 


5.801 ±0.2 


21 


6.0031 



Table 3.3: Frequency spacing between eigenmodes for the 7.32m long MH cavity pro- 
totype. The frequency values are expressed in MHz. The error is 0.2 MHz (sampling 
spacing and systematic effects). 
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Figure 3.26: Lowest order eigenmodes of the MH cavity with real deformed mirrors(both 
MH and Input). FEM eigenmodes calculation. 
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3.6.1.3 Tilt sensitivity 



Since the mesa beams are intended for use in actual interferometers it is important to 
study their ease of control. Some theoretical and experimental investigations have been 
carried out in this area |19( |20"1 [50j [62] . The alignment tolerances for a mesa beam arm 
cavity are expected to be ~ 3 times more stringent than those of the gaussian arm 
cavity in an advanced gravitational wave detector. 

In our 7.3 m cavity (just as in a long baseline GWID), extremely small tilts create 
significant modification of the mesa beam profile. 

The tilt precision of the optical lever is estimated to be 0.05 ^rad. 



Figure 3.27 shows the good agreement between our recorded profiles (thin lines) 
and FFT simulated data (thick lines). Note that these FFT simulated profiles were 
constructed using a two-mirror cavity (as opposed to the real three-mirror cavity). 



Mesa Beam deformation under MH mirror tilt 
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Figure 3.27: Comparison of FFT simulated (thick) and experimental profiles (thin) 
with same integrated power for various tilts. 
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Chapter 4 

Thermal noise in view of GW 
sensitivity improvements 



4.1 Introduction 

Thermal noise is one of fundamental noise sources in precision measurement, such as 
gravitational experiments [67] . It is expected that the sensitivity of the next generation 



interferometric gravitational wave detectors (see Fig. 1.3 ) will be limited by the thermal 
noise of its optical components (i.e. mirrors substrate and high-reflective dielectric 
coatings). Therefore, it is important to study the thermal noise in view of sensitivity 
improvement. In this chapter, the fundamental theorem and the estimation method of 
the thermal noise are introduced. Moreover, the evaluation of the thermal noise of the 
interferometric gravitational wave detectors are discussed. We will analyze two different 
method for thermal noise reduction: the first based on beam shaping whereas the second 
on coating optimization. In this thesis we calculate for the first time the coating noise 
reduction using mesa beam instead of standard Gaussian beam and we show that the 
mesa beam is a good candidate^ for application in advanced GW interferometers. 



4.2 Fluctuation Dissipation Theorem 

The fluctuation-dissipation theorem (FDT) is one of the most important theorems 
of non-equilibrium statistical mechanics. The FDT predicts the relationship between 

1 There are currently different studies on other beam shapes |68l I69j (and so other mirror 
shapes) to minimize the mirror thermal noise. We will comment briefly on this aspect later. 
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the spectrum of the thermal noise and the dissipation of lineaiQ dissipativd^ systems. 
Since the thermal fluctuations and the dissipation are governed by the same interaction 
between a system and the heat bath, there must be a relation between both phenom- 
ena. This theorem was established by Callen et al. [70 . In their fundamental work, 
they established a generally valid connection between the response function and the 
associated equilibrium quantum fluctuations, i.e., the quantum fluctuation-dissipation 
theorem. Another key development must be credited to Lars Onsager: via his re- 
gression hypothesis, he linked the relaxation of an observable in the presence of weak 
external perturbations to the decay of correlations between associated microscopic vari- 
ables |71t 172] . This all culminated in the relations commonly known as the Green-Kubo 
relations |73| [74] . This notion of "linear response" related to the fluctuation proper- 
ties of the corresponding variables constitutes the (response) fluctuation- dissipation 
theorem. 

For a macroscopic system, the fluctuation-dissipation theorem describes the be- 
havior of any generalized coordinate X of any system that is weakly coupled to a 
thermalized bath with many degrees of freedom. The interaction between the system 
and externals is represented by a fluctuating generalized force, F(t). The information 
of the dissipation of this system is included in the response of X to F. This response 
is derived from the equation of motion of the system. 

The information of the dissipation of the system is embedded in a generalized 
impedance function Z(u), and the FDT establishes a relation between the spectral 
density of the spontaneous fluctuations exhibited by the system in equilibrium and this 
generalized admittance (the inverse of the impedance) of the system. The impedance 
Z is defined as 

*<«) = ^ (4.1) 

where F and X are the Fourier components of the generalized force and velocity, 
respectively. If the system were completely conservative, then the impedance would 

1 The linearity states that the amplitude of the response is proportional to the amplitude of 
the driving force. 

2 A system is considered dissipative if, under an external solicitation, it absorbs energy 
through irreversible processes. 
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be perfectly imaginary. The bath prevents the system from being conservative: energy 
can be fed back and forth between the generalized coordinate X and the bath's many 
degrees of freedom. This energy coupling changes the impedance Z(uS) from pure 
imaginary to complex, and correspondingly, the resulting motions of X feed energy into 
the bath and vice versa. The FDT states that the spectral density (power spectrum) 
of the thermal fluctuation for the quantity X is 

S x (u) = 4 ^ft[Y{u)] (4.2) 

and that the spectral density of the generalized force Sf(uj) Qis given by 

S F (u) =4k b m[Z(io)} (4.3) 

The fluctuation dissipation theorem is derived from the assumption that the re- 
sponse of a system in thermodynamic equilibrium to a small external perturbation is 
the same as its response to a spontaneous fluctuation. Therefore it seems natural to 
look for a direct relation between the fluctuation properties of the thermodynamic sys- 
tem and its linear response properties to an external perturbation, which is coupled 
exactly to the generalized coordinate we are interested in. Levin |75] proposed a di- 
rect and very powerful method for the application of the FDT and derived a formula 
which requires no specific information about the dissipation mechanism, apart from 
the fact that it should be dissipation of mechanical energy and that the localization 
and the amount of dissipation are known. As a consequence, the Levin formula can be 
used to study a broad range of thermal noise induced fluctuations associated with all 
specific dissipations. We imagine to apply to the system an external oscillatory force 
F(t) = Foexp(iuJot) that derives the generalized momentum conjugate to X but does 
not derive the other generalized momenta of the system. In order to introduce this 
force into this equation of the motion, the new term, 

H int = ~F(t)X (4.4) 

is added to the Hamiltonian of the system. From the definition of the generalized 
impedance Z we have that X(t)) = z(u ) F(fy an d the average dissipated power within 

lr This is also known as generalized Nyquist noise formula 
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the system as a function of the driving frequency is given by 



(tWwo)) = (XF) = -Fffi[Y(u; )] 



(4.5) 



Substituting Eq. ( |4.5| ) into the expression (4.2) of the FDT, we get 



8k B T (W diss (co)} 



(4.6) 



where we have done the substitution ujq — > to since the same reasoning can be repeated 
for any driving frequency (in the range coupled to the bath) and the spectral density 
Sx has the claimed form anywhere in this range. This formula is particulary useful 
because for many systems, the calculation of the dissipated energy is simpler than that 
of the real part of the generalized admittance function. 

4.2.1 FDT: superposition principle 

In many practical situations it is necessary to evaluate the behavior of system in which 
the generalized coordinate we are looking at is a superposition of other generalized 
coordinates of the system. This can be the case of the surface of a mirror as seen by a 
laser beam profile, since the displacement sensed by the laser beam can be expressed as 
a spatially weighted superposition of the internal modes of the mirror interacting with 
the beam itself. Since the FDT is based on the linear response theory, it is obvious 
that for any linear combination of generalized coordinates we can repeat the reasonings 
of the previous section and obtain analogous formulas for the FDT. Lets define a new 
coordinate as a weighted superposition of certain system's generalized coordinates 



The generalized force F new that drives the momentum conjugated to X new is intro- 
duced in the Hamiltonian of the system, trough the following interaction term 




(4.7) 





and the generalized admittance of the system becomes 




(4.9) 
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where F and X new are the Fourier components of the generalized force and of the 
velocity of the new generalized coordinate respectively. The spectral density of the 
thermal fluctuation for the coordinate X new is given by 

Sx ne M = ^M[Y new (u)} (4.10) 

4.3 Mirror thermal noise 

The standard thermal noise is the phase noise caused by random motions of the reflect- 
ing faces of mirrors in a GW interferometer. A reflecting face can move either because 
it is displaced by its suspension systems or because it undergoes internal stresses. At 
finite temperatures both effects are possible. We address here the internal stresses. 
Consider a massive body at temperature T. If T > the atoms constituting the body 
are excited and have random motions around their equilibrium position. The fact that 
they are strongly coupled to neighboring atoms makes possible propagation of elas- 
tic waves of various types, reflecting on the faces and the onset of stationary waves. 
One can show that, for a finite body (like for instance a cylinder of silica) there is a 
discrete infinity of such stationary waves, each corresponding to a particular elastic 
normal mode. At thermal equilibrium, the state of the body can be represented by 
a linear superposition of all the modes, with random relative phases, and, due to the 
energy equipartition theorem, the same energy ksT {kg is the Boltzmann constant). 
The motion of atoms near a limiting surface of the body will slightly modify its shape, 
and if we consider the reflecting face of a mirror, a surface distortion is cause of phase 
change in the reflected beam; in other words, of a noise. Estimation of the resulting 
spectral density of phase noise in the probe beam is the internal thermal noise problem 
in massive mirrors. 

Modal expansion is a traditional method to calculate thermal noise. Although now 
superseded by the Levin approach, we find it useful to first discuss the general concept 
of this method and then discuss its application to the mirror. In a frequency domain, 
the equation of motion of a one-dimensional harmonic oscillator is written as 

-mu?X{uj) + mul[l + i</>(u)]X (u) = F{uj) (4.11) 
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where m is the mass of the oscillator, and cuq is the angular resonant frequency of the 
oscillator. We introduced the mechanical loss as the loss angle <f>(co), and mwo[l+i(^(w)] 
is called the complex spring constant. For a mechanical system with a spring constant, 
we describe the dissipation as the imaginary part of a complex spring constant. In the 
standard theory of inelasticity, Hookes law is extended to represent the internal friction 
process, leading to a description of the dissipation of energy in the oscillator in terms of 



its complex stiffness spring. As a matter of fact, Eq. (4.11 ) has merely phenomenological 
relevance to express the experimental appearance of a lag angle between the stress and 
the strain and the behavior with respect to the frequency. The mechanical impedance 
of the system is given by 



Z{uj) = — [<f>(u)u% + i(uJ 2 - ujI)} (4.12) 



and the application of the FDT (4.2) allows the calculation of the spectral density 



of the thermal fluctuations of the oscillator's position 

5xH = ^(^_4)2+°u; V 2 M (413) 
Before giving a short overview of the modal expansion approach to the calculation of 
the FDT, we have to summarize and fix the notation for the elastic theory of continuum 
media which will be used in the following Sections. 



4.4 Basic elasticity theory 

We label the position of a point (a tiny bit of solid) in an unstressed body, relative to 
some convenient origin, by its position vector x. Let a force be applied so the body 
deforms and the point moves from x to x + u(x); we call u the point's displacement 
vector. The strain tensor is defined as a second-rank tensor field given by the covari- 
ant derivative of the displacement vector which can be written in arbitrary system of 
coordinates using index notation and the symbol ';' for covariant derivative 

Sij = u i;j (4.14) 

In a Cartesian coordinate system the components of the gradient are always just 
partial derivatives. 
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The strain tensor S is a general, second-rank tensor. Therefore, its irreducible 
tensorial parts are its trace = Tr(S) = Si-i, which is called the deformed body's 
expansion; its symmetric, trace-free part S, which is called the body's shear; and its 
antisymmetric part R, which is called the body's rotation: 

The strain tensor can be reconstructed from these irreducible tensorial parts in the 
following manner 

S,j _',(->.'/,, • • />',, (4.15) 

Since the third term represent a rigid rotation which does not deform the solid and 
we are interested in the deformation itself, it is useful to rewrite the strain tensor as 

Sij = Eij + Rij = ^(ui-j + Uj-i) + ^(ui-j - uj-i) (4.16) 

The first part e is a symmetric tensor and represents elongation, compression, and 
shear. In many textbooks this part alone is called strain tensor. The displacement 
due to rotation will be not considered in the following since it is not related to the 
deformations of the elastic body. 

The forces acting within an elastic solid are measured by a second rank tensor, the 
stress tensor. Consider two small, contiguous regions in a solid. If we take a small 
element of oriented area cffi in the contact surface with its positive sense pointing from 
the first region toward the second, then the first region exerts a force dF (not necessarily 
normal to the surface) on the second through this area. The force the second region 
exerts on the first (through the area dXl) will, by Newton's third law, be equal and 
opposite to that force. The force and the area of contact are both vectors and there 
is a linear relationship between them. The two vectors therefore will be related by a 
second rank tensor, the stress tensor a: 

dF i = a ij dZj (4.17) 

and the elastic force per unit volume acting on a solid can be expressed as the 
divergence of the stress tensor 

f = 4 (4-18) 

The diagonal terms of the stress tensor define the normal stresses in the coordinate 
system in which is expressed. The off-diagonal terms define the shear stresses. 
Normal stresses act parallel to the unit vector orthogonal to the surface while shear 
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stresses act perpendicular to it. The conservation of angular momentum requires that 
the stress tensor be symmetric. Consequently, there are only 6 independent stress 
components to be determined. 

To understand the relationship between stress and strain in elastic media, we must 
generalize Hookes law for continuous, linear, elastic media. 

The most general linear equation relating two second rank tensors will involve a 
fourth rank tensor known as the elastic modulus tensor or elastic stiffness tensor, C. 
In tensor notation we have 

a ij = C ijkm £km (4 1Q) 

The inverse relation defines the fourth rank compliance tensor G 

£ij = Gij km a km (4.20) 

Now, a general fourth rank tensor in three dimensions has 3 4 = 81 independent 
components. However there are several symmetries that we can exploit. As the stress 
tensor is symmetric, and only the symmetric part of the strain tensor creates stress, 
C is symmetric in its first pair of indices and also in its second pair. Moreover, from 
the expansion of the elastic internal energy in terms of components of strain, it can be 
shown that the elastic tensor must be symmetric under the interchange of the first two 
indices with the second ones. There are therefore 21 independent components in C. 

Many substances, notably crystals, exhibit additional symmetries and this can re- 
duce the number of independent components considerably. The simplest, and in fact 
most common, case arises when the medium is isotropic. In other words, there are 
no preferred directions in the material. This occurs when the solid is poly crystalline 
or amorphous and completely disordered on a scale large compared with the atomic 
spacing, but small compared with the solid's inhomogeneity scale. 

If a body is isotropic, then its elastic properties must be describable by scalars. 
Now, the stress tensor a, being symmetric, must have just two irreducible tensorial 
parts, a scalar part and a trace- free symmetric part; and the parts of the strain that 
can produce this behavior are the scalar expansion and the trace-free, symmetric 
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shear X, but not the rotation. The only linear, coordinate-independent relationship 
between these tensorial quantities involving solely scalars is 

a tj = KQ 9ij + 2/iEij (4.21) 

Here K is called the bulk modulus and \x the shear modulus, and the factor 2 is 
included for purely historical reasons. Sometimes it is convenient to introduce the Lame 
coefficients A, \x which are related to K by the relation 

K = X+ 2 fi (4.22) 
3 

and write the elastic tensor in a generic coordinate system as 

C ljrs = \g lj g rs + ix{9 ir 9 3S + 3 J V r ) (4.23) 
and the relation between stress and strain tensor as 



a ij = 2fis ij + Xe k k g ij (4.24) 
In terms of the Young modulus Y and the Poisson coefficient v the linear elasticity 



relation (4.19) for an isotropic and homogeneous body can be written as 



a rs = (g mr g ns + g ms g nr + J^g™g m A e mn (4.25) 

1 + V \ 1 — Iv J 

If one has coordinate components, found from generalized coordinate tensor analy- 
sis, for some quantity, such as stress or strain, one needs to be able to translate those 
into the values measured in the experiment. The method for doing this is reviewed in 



Appendix A.l The measured value for a given vector component, unlike the coordinate 



component, is unique within a given reference frame. In differential geometry (tensor 



analysis), that measured value is called the "physical component" (App. |A.l ). It is 



important to recognize that physical components do not transform as true tensor com- 
ponents and one can not simply use physical components in tensor analysis as if they 
were. Typically, one starts with physical components as input to a problem. These are 
converted to coordinate components, and the appropriate tensor analysis carried out to 
get an answer in terms of coordinate components. One then converts these coordinate 
components into physical components as a last step, in order to compare with values 
measured with instruments in the real world. Since we will deal with only orthogonal 
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coordinate (Cartesian and Cylindrical) we can restrict our discussion to orthogonal 
coordinate systems in which the metric tensor is diagonal and 

ds 2 = gijdx l dx^ gu = hf, i not summed gij = i ^ j (4-26) 

The physical components of a tensor are given by the components of the tensor with 
respect to unit basis vectors. In the following sections we will deal with orthonormal 
cylindrical basis and the components denoted by r, (p, z are referred to this basis. More 
details about the relation between physical components and tensorial components is 
given in App. |A.l 



4.4.1 Thermal noise with mode expansion 

When the thermal noise of the internal modes of a mirror in interferometric gravitational 
wave detectors is calculated, the mirror is treated as a elastic cylinder. Cylindrical 
coordinates are employed with the origin at the center of the mirror and the z-axis along 
the cylindrical axis. The observable coordinate, X, which is the surface displacement 
averaged by the beam power distribution, is expressed as 

X = I u z (r)P(r)dS (4.27) 
Js 

where u z is the z component of the displacement vector u, and the weighting func- 
tion P is the normalized distribution of the beam power over the mirror surface. The 
external generalized force that drives the momentum conjugated to X is introduced in 
the equation of motion without dissipation for the elastic body in the following form 

d 2 u Y Y 

P^r--, ^Au-- r V(V • u) = F(t)P(r)z (4.28) 

H dt 2 2(1 + 2(1 + i/)(1 + 2i/) v ' w w y ' 

where p is the density and we have used the equation of motion of a linear and 
homogeneous elastic body (Eq. ( |A.14 | in App. A. 2 ). We can look for a solution of (4.28 ) 



in the form of a normal modes expansion, that is, we write the following superposition 
of basis functions w n 

u(r,t) =J2 w n(r)q n (t) (4.29) 
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The basis functions, w n , are the solution of the eigenvalue problem, for the free 
elastic body, written as 



-puj n w n 



Y 



-Aw„, + 



Y 



-V(V • w r 



2(l + i/)~"" ' 2(1 + i/)(1 + 2i/) ' v ' "' a ' (4 ' 30) 
where u; n and w n (r) correspond to the angular resonant frequency and the displace- 
ment of the n-th resonant mode of the system, respectively. 



The function, q n (t) in Eq.(4.29), represents the time development of the n-th mode. 



The equation of motion of q n is derived substituting the equation Eq. (4.291 for u in 



Eq. (4.28) and using the orthogonality property of the complete set of functions w n . 



The result is written in the form 



m n q n (t) + m n ujlq n (t) = F(t) 



(4.31) 



where the parameter m n , called effective mass of the n-th mode, is defined by the 
normalization choice for the functions w„ 



J w n , z {r)P(r)dS = 1 
m n = j p(r)\w n (r)\ 2 dV 



(4.32) 



where w n ^ z is the z component of w n . Therefore, the time evolution of the n-th 
mode is the same as that of a harmonic oscillator of mass m n and angular resonant 



frequency uj n with an external force F(t) acting on it. Using Eq. (4.29 ) and Eq. (4.32 )in 



Eq. (4.27) it follows 



X(t) 



E 



Qn{t) 



(4.33) 



showing that the observable coordinate X can be simply described as a superposition 
of the motions of the harmonic oscillators q n . Up to now we neglected dissipation in the 
system. To describe dissipation, the loss angles, <f> n (ui), are introduced in the frequency 



domain equation for the normal modes evolution. Equation (4.31) is rewritten as 



-m n LU 2 q n + ninOJ 2 [1 + i(f) n (u>)] • 



F. 



(4.34) 
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From Eq. (4.33) and the equations (4.34), the generalized admittance of the system 



is 



V ' F{u) F(u>) 



iuj 



E 



(4.35) 



-m n L0 2 + m n Un [1 + #n(^)] 
The power spectral density of X, Sx can now be derived applying the FDT, Eq. 



(4.10) and Eq. (4.35): 



S x (u) 



m r . 



Therefore the thermal motion of the system is expressed as the sum of the contribu- 



(4.36) 



tions given by the harmonic oscillators of the normal- mode expansion. Equation (4.36) 



allows to calculate the thermal noise from the angular resonant frequency uj n , the effec- 
tive mass m n and the loss angle (j) n (uj) of each mode. The angular resonant frequency 
and the displacement of the mode w n are obtained from the eigenvector problem, 
Eq. (4.30) . The effective mass m n is calculated from the w n found and equation 



(4.32). The loss angle is derived from the experiments but its measurement in a wide 



frequency range is commonly difficult; thus, it is usual to estimate it from the Q-value 
on the resonant frequencies for different modes, according to the following relation: 



Qr. 



1 



(4.37) 



4>n(u n ) 

There are two methods to solve the eigenvalue problem. The first one is the method 
proposed by Hutchinson [76J. This is a very accurate semi-analytical algorithm to 
simulate resonances of an isotropic elastic cylinder. The second method is the finite 
element method which is a numerical method. The thermal noise can be calculated 



using (4.36), where the summation has to be done over a number of modes that assures 



a good convergence of the series. Using the mode expansion formalism described above, 
the mirror thermal noise of interferometric gravitational waves detector has been calcu- 
lated by Gillespie-Raab [77] , Bondu-Vinet [78] under the assumption of homogeneously 
distributed loss inside the mirror bulk. 

Recently, a problem of the modal expansion method was clarified in [79]. The 
introduction of the mechanical loss after the system is divided into basis functions, is 
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not equivalent to the solution that is directly obtained from the equation of motion 
with the loss included since the beginning. The modal expansion fails, especially when 
the loss is inhomogeneously distributed in the system, (f>(r,u), because of a coupling 
between the internal modes. The LIGO mirrors, composed of several parts, mirror 
substrate, coating, magnet, each with its own loss mechanism, correspond to this case. 



4.5 Thermal noise: direct calculation 

At present the technique recognized to be the most appropriate to calculate the thermal 
noise is the "direct" method. It has the advantage that leaves aside the modal decom- 
position. There are many different ways of implementing this computation, Levins ap- 
proach [75], Nakagawas approach [80 j and Numata's numerical dynamic approach [81 . 
In the following we limit the discussion to Levin's method explained in Sec. |4.2| since 
this is the approach used in this thesis work. 

Consider the surface of the mirror invested by the laser beam which has a profile 
given by the weighting function -P(r). The read-out variable will be the generalized 
coordinate X expressed as a continuous combination of system coordinates in the form 



of Eq. (4.27). Using the FDT, Eq. (4.10), the calculus of the thermal noise reduces 



to the calculus of the real part of the admittance. Levin's method consist in the 
calculation of -KefY^^,] from the average mechanical energy dissipated by a particular 
driving force. Suppose to apply on the mirror surface an oscillatory continuous force 
F(r, t) that mimics the profile -P(r) of the laser beam: 

F(r, t) = F cos(wi)P(r) (4.38) 
This force corresponds to a generalized force 



F(t) = F cos(wi) 



(4.39) 



that drives only the momentum conjugated to X as is expressed in Eq. (4.4); then 



we can apply the equations (4.6) for the spectral density of the fluctuations of our 
observable coordinate X. 
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4.6 Mirror thermal noise contributions 

Mirror thermal noise (in the bulk and in the coating), can be divided into: 

• Brownian thermal noise, 

• thermoelastic noise, 

• thermorefractive noise. 

In this thesis, we will focus our attention on these kinds of mirror thermal noises 
and we will show the possible sensitivity improvements in using a mesa beam instead 
of a standard Gaussian beam. 

4.6.1 Brownian thermal noise 

Internal friction in solids was identified by Kimball and Lovell [52], who described 
it as a phase shift between stress and strain. Brownian thermal noise due to mir- 
rors, can be interpreted as a fluctuation of the mirror surface position coming from 
the mirror modes thermal excitation or, using the FDT, as fluctuations induced by a 
structural damping. In the real lattice, impurities, dislocations and imperfections play 
the role of dissipation sources in the phonon dynamics. Irreversible processes can also 
be attributed to relaxation mechanism in asymmetric potential- well models [83} 184"] . 
Thermal noise due to homogeneously distributed damping processes, such as Brownian 
bulk noise in GW detector mirrors, can be estimated using complex valued elastic co- 
efficients whose imaginary part is related to the dissipation mechanism. For example, 
the structural damping can be parameterized by a complex materials Young modulus, 
were the imaginary part is related to the energy dissipation of the system 

Y\ — >Y(l + i<j>) (4.40) 

Measurements [85] have shown that 4> is dependent on the frequency /, but its 
dependence is sufficiently slow that it can be neglected in the frequency region of 
interest. To calculate the power spectral density associated to Brownian bulk motions, 



Levin expressed (Wdiss) m Eq. (4.6) as 



(Wdiss) = 2u;(f> s (U s ) (4.41) 
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where tt s is the loss angle of the substrate, (..) denotes the time average over the 
period of the oscillatory pressure previously defined and U s is the elastic energy stored 
in the substrate of the test mass. This energy is derived by using elasticity theory, and 
will depend on which monitoring beam profile you use as well as the mirrors dimensions. 
Levin supposed that the mirror is an infinite half-space and that the elastic deformation 
follows quasistatically the oscillatory pressure (a reasonable hypothesis because the first 
resonance mode of the mirror is far higher than the region of interest). In the case when 
the beam has a gaussian profile and the mirror is approximated by a semi infinite half- 
space, the displacement spectral density will be [73] 



The calculation was also made by Bondu et al. [86 , and reviewed by Liu and 
Thorne [57] in the case of a finite size mirror. 

The correction to the noise spectral density is an infinite series of Bessel function 
with coefficients depending on the material properties and on mirror and beam-spot 
size dimensions. If compared with the result in the approximation of infinite mirror 
mass, this calculation gives a correction to the thermal noise that, e.g. in the case of 
the new Advanced LIGO mirrors, starts to be greater than 30% for beam spot sizes 
greater than 6 cm. 

Considering that there is the project to use flat beams, this correction becomes 
crucial in estimating the Brownian thermal noise correctly. 

4.6.2 Thermoelastic noise 

Thermoelastic noise is intended as a noise that comes from the coupling of thermal 
fluctuations with displacement fluctuations thanks to a non-null coefficient of thermal 
expansion. The oscillating squeeze and stretch of the substrate material causes an 
oscillating, inhomogeneous temperature distribution: heat flows down the temperature 
gradient in such a way that it converts oscillation energy into additional heat. 

The inhomogeneity of the deformations causes the temperature perturbation 8T 
to be inhomogeneous, and that inhomogeneity produces a heat flux q = —kSJST. 
Whenever an amount of heat flows from a region of temperature T to one of slightly 



S 



■B,s 



Ak B Tl-v 2 



(4.42) 
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lower temperature T — dT , there is an increase of entropy and the resulting rate of 
entropy increase per unit volume is 

dS -q-V^T _ k(V5T) 2 
dVdt T 2 T 2 { ' 

This entropy increase entails a creation of new thermal energy at a rate per unit 
volume dEth/{dVdt) = TdS/(dVdt). Since, for our thought experiment with tempo- 
rally oscillating applied stress, this new thermal energy must come from the oscillating 
elastic energy, the rate of dissipation of elastic energy must be 

W dlss = J ^^dV (AAA) 

To calculate this noise for the mirror it is necessary to solve a system of two coupled 
equations, the first one is the elasto-dynamic equation for the displacement u(r,i) 
including a stress term coming from the temperature inhomogeneity, and the second 
one is the thermal conductivity equation for the temperature perturbation ST(r,t) 
including an heat-source term coming from the non uniform expansion. As the time 
required for sound to travel across the mirror is usually smaller than the oscillatory 
period, it is often used a quasistatic approximation in the equation of elasticity for a 
field of deformation u 

V(V • u) + (1 - 2z/)Au = 2a(l + u)V5T (4.45) 
The thermal conductivity equation for the temperature inhomogeneity 5T 

dST Kj aYT 0( V • u) , . 

A<5T = (4.46) 

dt C C(l-2v) dt K ' 

where Y, u, a, k and C are Young modulus, Poisson ratio, the coefficient of linear 

thermal expansion, the thermal conductivity and the specific heat per unit volume 

at constant volume of the mirror. Braginsky et al. [88] solved the problem using the 



fluctuation-dissipation theorem in the form (4.10), approximating the mirror with an 
infinite half-space and considering that the oscillatory period is far higher with re- 
spect the typical time scale of the diffusive heat flow. This is the so called adiabatic 



limit in which the second term on the left hand side of (4.46) is neglected. In this 



approximation, 5T is simply proportional to the elastic expansion 0. 
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Liu and Thorne |87J solved the same problem with the same approximations ap- 



plying the FDT in Levins form, Eq. (4.6), and using (4.44). They found that the 



displacement power spectrum for half infinite mirror is 



~,TE,s 



16 2/-, , \2 kuT K 1 

—a (1 + v) 3 



(4.47) 



They also made the calculations for a finite cylindrical test mass finding the cor- 
rection factor w.r.t. the infinite test mass solution. Depending on the mirror/beam 
geometry and the mirror's material this correction can be several 10%. 

Cerdonio et al. [89] extended the previous analysis releasing the assumption of the 
adiabatic limit, i.e. performing the analysis that is valid also at low temperatures 



and/or small beam spot size. Equation (4.47) has to be multiplied by a correction 
factor, which depend on the ratio between the beam radius and the thermal diffusivity 
length (that is frequency dependent). 



4.7 Coating thermal noise 

It has been pointed out [75] that it is critically important how the losses are distributed 
inside the test masses. Losses far from the beam spot contribute less to the total 
thermal noise, whereas losses near the spot, for example in the dielectric coating directly 
reflecting the beam, contribute more. 

The coatings on the mirrors play the crucial role in reflecting the laser used for 
position sensing off of the test mass. They must be highly reflective and able to handle 
high optical power, to reduce shot noise, and have low levels of intrinsic noise. Primary 
among intrinsic noise are thermally driven motions of the coating face and the optical 
path length. 

To obtain the required high reflectivity, multi-layers, dielectric coatings are used. 
Such coatings consists of alternating ion-beam deposited layers of optical thickness 
equal to A/4 (at 1.064// m) of two dielectric materials with differing refractive indexes 
(i.e., Ta<iO^I S1O2 for initial LIGO). It was found [90] that the mechanical loss comes 
from internal dissipation of the coating materials, and it is the high index Tantalum- 
pentoxide that is primarily responsible for the thermal noise in first generation inter- 
ferometers. 
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Another source of thermal noise is thermally driven fluctuations in the optical path 
length of the coating. This can occur from thermal fluctuations in length, i.e. thermoe- 
lastic noise |91( 192] . or from fluctuations in the refraction index, i.e thermorefr active 
noise [93]. These two mechanisms are driven by the same thermal fluctuations, there- 
fore will add coherently. This collective noise source is referred to as thermo-optic noise. 
The level of thermo-optic noise in next generation gravitational-wave interferometers 
depends on the thermal, elastic, and optic properties of the coating materials, primarily 
on a(= dL/dT) and (3{= dn/dT). Values in literature j9T| [92]for ion beam deposited 
coatings indicate that thermo-optic noise could contribute to Advanced LIGO sensitiv- 
ity and that the high index coating material may be the primary problem and lots of 
effort by the optics research groups is dedicated to measure these parameters. 

4.7.0.1 Brownian coating thermal noise 

Coating thermal noise due to internal losses is expected to be the dominant contribution 
to the thermal noise for mirrors with a SiOijTaiO^ coating on a fused silica substrate. 
Nakagawa et al. [80 calculated the dissipation induced by an inhomogeneous loss dis- 
tribution on a half-infinite mirror due to the fact that loss angle in the bulk material 
<j) s is different from loss angle in the coating <j) c but neglecting the differences in the 
elastic properties of the substrate and the coating 

gB = wi^'tMM' (4 . 48) 

y7T UjY s W yJlT (1 — V) W y 7T W 

Harry et al. [93], using Levins method, generalized this result for a non isotropic 
coating loss, assuming that the loss angle 0ii associated with energy stored in strains 
parallel to the plane of the coating, is different from the loss angle <j>± associated with 
the energy stored in strains perpendicular to the surface. For a thin coating, provided 
that v c ~ Y c = Yj| ~ Y± and u c = v\\ ~ u±, they found 



lB>c = 4k B Tl-v 2 1 d y c (l + u s ){l - 2v s ) 2 + Y s v c (l + v c )(l - 2v s 
' x ojYw Y s (1 - !/2)(l - „,) 



Y s (l + u c )(l - 2u c ) - Y c v c {\ + v s ){\ - 2v s ) \ 

+ nd-,c)(i-^) ^ (449) 
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Equation (4.49 ) is valid provided that the coating loss occurs in the coating materials 
themselves and not for rubbing in the bulk/coating or layer /layer and this has been 
experimentally investigated by Penn et al. |90( [95] . 

4.7.0.2 Coating thermoelastic noise 

Temperature fluctuations (that couple with displacement thanks to a can be originated 
either by intrinsic thermodynamical fluctuation, or by laser-photon absorption. We will 
briefly describe the first mechanism and refer to Rao [96] for the description of the latter 
one since that is a negligible effect in LIGO and Ad-LIGO mirrors. 

The diffusive heat characteristic length ry of the substrate and coating (of the order 
of mm) is far larger than the coating thickness (a few fi m). Because diffusive heat 
flow in the longitudinal direction is not negligible, heat flow in the direction normal to 
the coating cannot be treated adiabatically. The problem was extensively treated by 
Braginsky and Vyatchanin [£JT] and by Fejer et al. [92 in the case of Gaussian beam. In 
the frequency band f < t^ 1 , where r c is the thermal diffusion time across the coating 
(t^ 1 ~ 6 kHz for TCI2O5/ ' S1O2 multi-layer coating), coating thermoelastic noise can be 
estimated by the equation 



where A 2 is a dimensionless combination of material constants that vanishes when 
the film and substrate are identical 



A 2 »' C» a < 



2a,C r (l 



rr^ + (1 - 2 ^% 



1 (4.51) 



4.7.0.3 Coating thermorefractive noise 



The thermo-refractive noise is generated by temperature fluctuations that couple with 
phase fluctuations of the reflected laser beam (and therefore with measured displace- 
ment), thanks to the non- vanishing coefficient /3 = where n is the refractive index. 
Considering the frequency range (that is typical for the ground interferometer) around 
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Thermal noise 


S x oc w~ n 


Substrate Brownian 


n = 1 


Substrate Thermoelastic 


n = 3 


Coating (all kinds) 


n = 2 



Table 4.1: Dependence of the thermal noise spectral densities on the Gaussian beam 
radius. 

100 Hz and an half infinite mirror with a Gaussian beam, Braginsky et al. [93] calculated 
the equivalent displacement noise spectral density 



_ 2V2 k B //3 e/ /AT x 2 

7T yJCKUJ V w 



where /3 e / / = ~^~irz^y f° r a multi-layer coating which consists of alternating sequences 
of quarter-wavelength dielectric layers having refractive indices n\ and njfj] 



4.8 Thermal noise: Gaussian VS mesa beam 



Coating thermal noise (4.49) due to internal losses is expected to be the dominant con- 
tribution to the thermal noise for mirrors with a SiO^lTa^O^ coating on a fused silica 
substrate, whereas the thermoelastic noise of the substrate is the dominant contribu- 
tion for sapphire mirrors at room temperature. The local surface fluctuations produced 
by thermal noise are averaged by the intensity distribution of the laser beam spot over 
the mirror surface. Reading the entire mirror surface with uniform sensitivity would 
minimize the thermal noise. The standard design of interferometers uses light beams 
with a Gaussian distribution of power, which are eigenfunctions of cavities with spher- 



ical mirrors, a well-developed and understood technology. Table 4.1 summarizes the 
geometrical dependence of each kind of mirror thermal noise on the Gaussian beam 
radius w for half infinite test mass. 

The larger is the beam radius, the better is the averaging of the fluctuations and 
thus lower will be the noise. However the beam size is constrained by the allowable 

1 This formula has been revisited many times during the last few years. Originally was 
proposed in this form in |93j ; subsequently was modified in |91j . and recently an independent 
analysis |97j confirmed the original formula. 
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diffraction losses requirements, which cannot exceed a few ppm. Taking into account 
the diffraction loss constraints, a Gaussian beam, effectively averages out the thermal 
fluctuations only over a few percent of the mirror surface. A significant reduction 
in all kinds of mirrors thermal noises can be achieved by using modified optics that 
reshape the beam from a conventional Gaussian profile into the "mesa-beam" profile 
described in Sec. |3.4| A large-radius, flat-topped beam with steeply dropping edges 
(necessary to satisfy the diffraction loss constraint) will lead to a better sampling of 
the fluctuating surface, lower noise in the determination of the mirror surface position 
and better sensitivity for GW detectors. The calculation of substrate thermo-elastic 
noise reduction using Mesa beam has been done in |50] for sapphire test mass. More 
recently Vinet [52] calculated the substrate Brownian thermal noise reduction using 
Mesa beam for Virgo mirror size. 

Calculation of the coating thermal noise, which is expected to be the most significant 
contribution to the thermal noise budget for the test masses of the next generation 
of GW interferometers, has never been published for non Gaussian beams and finite 
cylindrical test masses. In this thesis we present a comparative study of the various 
sources of thermal noise in different mirror and beam configurations, considering both 
Gaussian and Mesa beam profiles, addressing the problem of thermal noise reduction, 
through mirror size aspect-ratio and beam size optimization. (Some of these results 
have been already presented in [5T] and [53 ]). 

4.8.1 Elastic solution for a cylindrical test mass 



We used the Levin approach to the Fluctuation Dissipation Theorem (4.6) in order 
to calculate the power spectral density of the test mass displacement. The generalized 
coordinate X is given by the average of the normal displacement of the test mass surface, 



weighted by the beam spot power distribution as introduced in Eq. (4.27). In this way 
the thermal noise evaluation reduces to the calculation of (Wdiss) i n accordance to the 
specific mechanism and localization of dissipation, i.e. Brownian or thermoelastic and 
substrate or coating. 

The main ingredient for our calculations is the model proposed by Bondu, Hello 
and Vinet [86] and corrected by Liu and Thorne [57] of the approximate solution of 
the elasticity equations for a cylindrical test mass subjected to the oscillatory pressure 
with the same spatial profile as the beam power distribution (assumed cylindrically 
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symmetric). Consider a finite sized, cylindrical test mass with radius a and thickness 
H and with a cylindrically symmetric light spot centered on the cylinder's circular face, 
which applies an oscillating pressure 

P(r) = F /(r) cos(wi) (4.53) 

where Fq is a constant force amplitude and /(r) beam power distribution normalized 
over the mirror surface 

2vr / /(r)rdr = 1 (4.54) 
Jo 

The fundamental approximation underneath these calculation is the quasi-static 
approximation for the calculation of the displacement fields according to the oscillatory 
pressure which is a good approximation for oscillatory period larger than the time 
required for sound to travel across the test mass and far away from the resonances of 
the test mass. 



The static stress balance equations in cylindrical coordinates derived in Appendix A. 2 
have to be solved together with the following boundary conditions: 

• No shear on the cylinder's surfaces, i.e. 

ffrz(r = a> z) = 0, cr rz (r, z = 0) = 0, o rz {r, z = H) = (4.55) 

• Pressure of the beam on the front face and no normal stress at the back surface 



cr zz {r, z = 0) = -P(r), a zz {r, z = H)=0, (4.56) 
• No radial stress on the cylindrical edge 

a rr (r = a, z) = (4-57) 

The displacement vector satisfying the mentioned boundary conditions Q is given 



by 



1 The boundary condition (4.57) is not fulfilled exactly, in fact, as discussed in [HS], the Co 
and c\ terms in the displacement arc a correction to the leading-order displacement, designed 
to improve the satisfaction of the a rr (r = a, z) = boundary condition as far as a rr (r — a, z) 
calculated without these corrections can be approximated by a linear function. 
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u r (r, z, t) = Fq cos(ujt) 



^2 A m {z)Ji{k m r) + n ^ 2fi n ^ (c r + c x r z) 



m=l 



+ 



2 / u(3A + 2/x) 



Xp r 



(-1) 



(4.58) 



ue(r,z,t) = 



(4.59) 



u z (r, z, t) = Fq cos(wi) 



A / c z 2 
^B m (z)J (k m r) - M3A + 2/i) + ^ 



A + 2,u 2 

cir + 



Ap r 2 



4/n(3A + 2fi) 

(A + ^)p 
/x(3A + 2/x) 



2F 



(4.60) 



where Jo and Ji are the Bessel function of order zero and order one respectively, 
Cm is the m— th zero of the Bessel function J\(x), k m = Cm/a and the other coefficients 
are 



Po 



B m (z) 



7TGT 



Pn 



" 2 -'o(C/»J ./() 



f(r)J (k m r)rdr 



a 2 s-^ Jo((m)p m 
c o = 6 772 2^ — 72 > 

Sm 



Cl = - 



m=l 



2co 



2 A + 2// 



A m (z) = 7m e- fc ™ z + 5 m e fem2 + ( a m e- fcm2 + (3 m e kmZ ^j 



k m z \ + fi 
A + 3/U 



(a m e kmZ - /3 m e fc ™ 2 ) + 



A + 3fi 
2(A + 2/x) 



/3 m (5^ 



[2(X + 2ft) 

and a m , P m ,^ m ,S m are constants given by 



(4.61) 
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Qm = e~ 2kmH (4.62) 

Pm 

(A + 2fj) 1 - Q m + 2k m HQ m 



Oim 
Pm 
7m 



k mf i(\ + m) (i - Qm) 2 - ^ m mQ m 

Pm 

(A + 2fi)Q m l-Q m + 2k m H 
k m ^{\ + fi) (1 - Q m ) 2 - 4k m H 2 Q m 

Pm [2k 2 a H 2 {\ + n) + 2fik m H]Q m + fijl-Qm) 
2k m fi(\ + fi) (1 - Q m y - 4k 2 m H 2 Q m 



r PmQm 2k 2 n H 2 (X + fi) -2fik m H - - Q m ) 

m ~ 2k m /j l (\ + //) (l-Q m ) 2 -4fc^ Qm ' j 

In the axi-symmetric case, the non-zero (physical) components of the strain tensor 
in cylindrical coordinates are given by 

du r u r du z 1 fdu r du z \ 

Or r oz 2 \ Oz Or J 

and, for an homogeneous and isotropic body, the stress components aij are related 

to the strain by the Lame coefficients A, /j, 



c rr — AG) -|- 2fis rr 
aee = A0 + 2/ieee 
Ozz = A0 + 2fie zz 

a rz = 2[ie rz (4.65) 
where G = e rr + egg + e zz is the expansion. 

4.8.2 Thermal noise calculations for finite size test masses 

In order to manipulate these quite complicated expressions, we developed a Mathematica® 
Notebook called Thermal Noise Notebook (TNN)Qwith the intent of providing a user- 
friendly tool which can be used by GW interferometer researchers to evaluate quickly 
the expected mirror thermal noise. We chose to implement the calculation using this 
program because of its great versatility for numerical and analytical calculation and 
the possibility of changing very easily all the parameters involved in the calculations, 
1 Available at http : //www . ligo . caltech . edu/~jagresti/| 
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beam shape, mirror's aspect ratio, mechanical and thermal properties of the material of 
the substrate and coating. We now give an overview of the calculation for the different 
types of thermal noise. 

4.8.2.1 Substrate Brownian thermal noise 



The conventional thermal noise of the substrate is given by Eq. (4.6) with the time 
averaged dissipation (Wfc 8a ) given by 

(W diss ) = 2u;<t> s (U s ) (4.66) 

where (f) s is the loss angle of the substrate,(..) denotes the time average over the oscilla- 
tory period and U s is the elastic energy stored in the test mass that can be determined 
integrating, over the test mass volume, the elastic energy density given by 

g(r, z) = -ei^ = - (A9 2 + 2fi(e 2 rr + e 2 ee + e 2 zz + 2e 2 TZ )) (4.67) 

U s = I g(r,z)dV (4.68) 
Jv 

e%j and Gij are the component of the strain and stress tensor respectively calculated 



from the displacement vector using equations (4.64) and from the constitutive relations 
(elastic moduli tensor) for an homogeneous and isotropic material, Eq. (4.65). These 
cumbersome calculations are performed very efficiently by Mathematica® and will not 
be reported here. 

4.8.2.2 Substrate thermo-elastic noise 



As reported in Sec. 4.6.2 the thermo-elastic dissipation is given by 

\2 



(W diss ) = I Ks{V J TsY dV (4.69) 



Where k s is the substrate thermal conductivity and 5T S is the substrate temperature 
perturbation induced by the elastic deformation due to the oscillatory pressure and is 
given by 
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Where a s is the linear thermal expansion coefficient, Y s and v s are the Young 
modulus and Poisson ratio respectively, C s is the specific heat per unit volume of the 
substrate. 



The equation (4.70) follows from the adiabatic approximation of the general ther- 



mal conductivity equation (4.46). The adiabatic approximation, already discussed in 



Sec. |4.6.2 consist in the following: if the time scale for diffusive heat flow is much longer 
than the pressure oscillating period, we can approximate the oscillations of stress, strain 
and temperature as adiabatic, neglecting the heat flow term in the thermal conductivity 



equation. Using equations (4.64) for the calculation of the expansion 0, substituting 



in (4.70) we can calculate the thermo-elastic dissipated energy in (4.44) and then the 



spectral density of the displacement noise given by Eq. (4.6). 



4.8.2.3 Coating thermal noise 

The geometry we consider for the reflective surface of the mirrors consists of a thin 
film of thickness on a substrate whose thermo-mechanical properties are different from 
those of the film. To simplify the analysis, we assumed that the multi-layer coating can 
be approximated as a uniform layer with appropriately averaged properties. 

4.8.2.4 Coating Brownian thermal noise 

In this case the averaged energy dissipated by the intrinsic losses in the coating is given 



by an analogous formula of (4.66) 



(W dlss ) = 2u>MU c ) (4.71) 

Where cp c is the loss angle of the coating and U c is the portion of elastic energy 
stored in the coating (in this calculation we assume an isotropic and homogeneous 
coating with averaged elastic coefficient: Young modulus Y c and Poisson ratio v c ). In 
the thin film approximation d <C H, we assume that the energy stored in the coating 
is given by U c = d5U c , where d is the thickness of the coating and 5U C is the energy 
density stored at the surface, integrated over the surface, which is expressed in terms 
of the coating stress and strain as 



2 £ ij a ij 



(4.72) 
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In the thin film approximation, the coating is approximated as a thin layer in 
which the stress and strain do not vary considerably as a function of depth within the 
coating and we shall approximate them as being constant. Following [94 , the stresses 
and strains in the coating can be calculated in terms of the stresses and strains at 
the surface of the substrate because of the boundary condition between coating and 
substrate: the coating must have the same tangential strains as the surface of the 
substrate and the coating experiences the same perpendicular pressure as the surface 
of the substrate. Since the only exerted force is normal to the plane we must have no 
shear stress on the coating, a c rz = 0,(a^ z (r, z = 0) = is a boundary condition for the 
elastic problem of the substrate) . The coating stress and strain can thus be calculated 
by the following equations 

e c rr (r) = e s rr (r, z = 0), e c ee (r) = e s ee {r, z = 0), e c zz (r) = e s zz {r, z = 0) 

<r% k = A c e c + 2fi c e c kk with k = r,0,z and 9 C = e c „. + e% e + e c zz (4.73) 

In this way we can calculate all the fields necessary for the computation of the 
elastic energy stored in the coating, U c , using the expressions already found for the 
substrate. 

4.8.2.5 Coating thermo-elastic noise 

In the thermo-elastic problem of the coating is important to note that the coating 
thickness, the diffusive heat transfer length and the beam radius satisfy the following 
relation 



d«.r T <.w, with r T = J ^ (4.74) 

This relation justify the approximation of the multi-layer film as a uniform film with 
averaged properties and when computing the oscillating temperature distribution we 
can consider the temperature variation as adiabatic in the transversal direction and that 
only the thermal diffusion orthogonal to the surface of the mirror need to be considered. 



The two thermo-elastic coupled equations (4.45) and (4.46) are solved perturbatively 



at the first order in the linear thermal expansion coefficient a. We first consider the 
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quasi-static stress-balance equation at the zeroth order in a which has been solved in 



the preceding section (Eqs. 4.73). Then, we solved the thermal conductivity equation 



Eq. (4.46) for a one- dimensional heat flow and using the zeroth order elastic fields as 



the source term. 



d d 2 \ 



YaT d9(z = 0) 
(1 - 2v)C dt 



P 



P 



(4.75) 



where Q(z = 0) is the expansion at the mirror surface associated with the zeroth- 
order elastic fields calculated for the previous sections and (5 = s,c indicates quantities 
evaluated in the substrate and the coating respectively. For a multi-layer coating this 
equation determines an averaged temperature field and the coating quantities are aver- 
aged following [92] . If d\ and d 2 are the thickness of the two materials composing the 
coating (d\ + d 2 = d) , we have 



P0 



avg 



di +d : 



dx + di 



/S.-y 



t\' r . 



-1 



(CB) 



avg 



(4.76) 



ft 1 '"C V" iwutjt — t pi 

^avg ^avg 

Assuming a time dependence of the form e lu)t for the oscillatory thermal and elastic 



fields, equations (4.75) can be cast in this form 

d 2 



iu> — K 



ST, 



P 



iuB 



P 



s, c 



(4.77) 



with the boundary conditions of zero heat flux at the surfaces of the test mass and 
continuity of temperature and heat flux at the boundary between coating and substrate 



d5T c 



dz 



dST s 



2 = 



dz 



0, 5T C = 8T S 



z=H 



z=d 



The general solution of (4.77) is given by 



05T C 
dz 



d5T„ 



dz 



z=d 



5Ta 



-Bp + Clpe^ z + C2 P e~^ z , 7/9 = (l + i) 



2K, 



P 



(4.78) 



(4.79) 



The boundary conditions (4.78) determine the four arbitrary constants Clg,C2g 
and the solutions for the temperature variation induced in the coating and in the 
substrate are 
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5T C = -B c + e T^( e 2 7ad _ e ^H) Ks ( Bc _ Bg)ls [ e7c2 + e - 7c ^j . (4 _ g0) 
^ h (k c1c - k s1s ) + e 2d ^ + ^\-K clc + k s + 7s )+ 



ST S = -B s + e ^ d (e 2 ^ d - 1)k c (B c - B s ) lc [e^ z + e 2H ^ e~ • (4.81) 
^ h (k c7c - k s7s ) + e 2d ^ + ^(-K clc + k s + 7s )+ 

+ e 2 ^ d ( Kc7c + k s7s ) - e 2d ^ +2H ^( Kc7c + k s + 7s )j 

The dissipated power for the coating thermo-elastic noise is given by 

k s ( d5T s \ 2 f k c ( dST c x 



w *"*J v .nnr) dV ' + J v f{^f) iV < (482) 

where we are neglecting the dissipation due to the radial heat flow [^] and the integrals 



are extended over the substrate and coating volume respectively. Averaging Eq. (4.82) 



over the oscillatory period, and inserting the result in Eq. (4.6) we have the spectral 



density of displacement noise due to coating thermo-elastic fluctuations. 
4.8.2.6 Coating thermo-refractive noise 

Thermo dynamical fluctuations of temperature in mirrors of gravitational wave antennae 
may be transformed into additional noise not only through thermal expansion coefficient 
but also through temperature dependence of refraction index. In this case we cannot use 
the original Levin's method^] to calculate this noise because it is not associated with 

lr The correction is completely negligible for the coating contribution, and is around 3% for 
the substrate dissipated power. 

2 Recently Levin proposed a variation of its original direct method, to calculate the thermo- 
refractive noise [arXiv:0710.2710]. 
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mechanical energy dissipation. Following [93], the spectral density of the equivalent 
displacement noise can be written as 

SxH = \ 2 (3 2 eff S T (u) (4.83) 

where A is the laser wavelength, St is the spectral density of temperature fluctuation 
in the coating and /3 e // for a multi-layer quarter wavelength optical coating is given by 

Peff = i^^f (4-84) 

The spectral density of the temperature fluctuation can be calculated using the 
Langevin approach explained in [93 . For semi infinite mirror and for cylindrically 
symmetric beam profile fir), it is easy to derive the following equation 

4k B T 2 K [°° f°° q ± dq ± 2q 2 1 ~ 2 

St{UJ) = —^L dqZ J (2vr) 2 WT^H^ l/(9i)l 

Q 2 = Q± + ll = Ql + Q 2 y + Qz (4-85) 

where f(q±) = 2tt rdrf(r)Jo(q±r) is the Hankel transform of the normalized 
power distribution over the mirror surface. 

In this section we want to give just an estimation of the coating thermo-refractive 
noise reduction using a fiat beam instead of a Gaussian beam and for this purpose we 
will approximate the real mesa beam as a perfect Flat Top beam 

r < b 

(4.86) 

r > b 

For Gaussian and Flat Top beam the Hankel transform can be analytically per- 



formed but the integral in (4.85) must be numerically evaluated (unless other approx- 
imations are accepted, as done in REF). For this comparison we chose a value of 
b = 4:Wq, (wo = 2.6cm) which correspond to the "standard" radius of the integration 
disc for Mesa Beam and compare this ideally fiat beam with the Advanced LIGO Gaus- 
sian baseline design w ~ 6cm. The displacement noise is reduced by a factor 1.7 in the 
case of a Flat Top beam. 



qFT 



~ 1.7 (4.87) 

f=100Hz 
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4.8.3 Results and discussion 

This is a scheme of the calculations performed with TNN for the analysis of the expected 
mirror thermal noise in Ad-LIGO interferometer: 

• we fixed the mirror mass at 40 Kg, constrained by the Advanced LIGO suspension 
system design 



We fixed the diffraction loss constraint at lppm(10~ 6 ) for both the Gaussian 
and the Mesa beam, calculating the diffraction losses with the so called clipping 
approximation; in this approximation the losses are computed by the amount of 
light that falls outside the mirror and the beam profile is assumed to retain its 
shape even though the diffraction from the edge of the mirror. 

We changed the mirror radius from 12 cm to 21 cm (with 1 cm per step) and at 
the same time we increased the Gaussian beam radius w and the Mesa beam 
integration disc radius b to satisfy the lppm constraint for diffraction losses. The 



corresponding values are shown in Fig. 4.1 The thickness H of the mirror is 



reduced correspondingly to satisfy the total mass constraint. In this way all the 
geometric parameters in the problem are functions of the mirror radius a. 



4.8.3.1 Fused silica substrate 



Figures 4.2(a), 4.3(a) |4.4(a) 4.5(a) show the displacement noise of all the analyzed 
thermal noise contributions, for Gaussian and Mesa beam in the case of Fused Silica 
substrate. The dominant noise, the coating Brownian thermal noise, undergoes a reduc- 
tion of a factor 1.7 for a mirror radius of 18 cm. The substrate thermal noise is reduced 
by a factor 1.55 , whereas the coating thermo-elastic and the substrate thermo-elastic 
are reduced by factor 1.7 and 1.9 respectively. 
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Figure 4.1: Mesa beam and Gaussian parameters for Ad-LIGO like cavity. 
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(a) Fused Silica (FS) substrate. 
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Figure 4.2: Substrate Brownian (SB) thermal noise. 
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(b) Sapphire (S) substrate. 
Figure 4.3: Substrate thermoelastic (ST) thermal noise. 
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(b) Sapphire (S) substrate. 
Figure 4.4: Coating Brownian (CB) thermal noise. 
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(b) Sapphire (S) substrate. 
Figure 4.5: Coating thermoelastic (CT) thermal noise. 
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4.8 Thermal noise: Gaussian VS mesa beam 



It is interesting to note that all the thermal noise contributions present a minimum 
in the finite cylindrical model. The minimum of the sum of all contributions represents 
the best choice for the mirror and beam dimensions. The descent part of the noise 
curves reflects the basic idea that increasing the beam radius the noise will get lower 
and the rising part of the curves can be explained heuristically by the fact that all the 
noises contributions are related somehow to the amplitude of the elastic deformation of 
the test mass under a surface pressure and this effect is larger in gong-shaped mirrors 
than in bar-shaped ones. 

From the incoherent sum of the four contributions we have that minimum thermal 
noise occurs for a mirror radius of about 18 cm for Gaussian beam and for about 19 
cm for Mesa beam. The corresponding mirror aspect ratios (2a/H ) are 2 for Gaussian 
beam and 2.4 for Mesa beam. The aspect ratio chosen for Advanced LIGO baseline 
is 1.7. For this geometry the gain in sensitivity is about a factor 1.7 switching from 



Gaussian to Mesa beam at the minima of thermal noise (see Fig. 4.7). 



Fig. 4.6 shows the expected sensitivity of Advanced LIGO interferometer in the 
Gaussian or Mesa beam configuration calculated with 5enc/i[[] The estimated range 
for NS-NS binary systems increases from 177 Mpc for Gaussian beam to 228Mpc for 
Mesa beam. This is a remarkable factor if we consider that we didn't optimize the 
other interferometer's parameters to take full advantage of the reduced mirror thermal 
noise floor. In this evaluation the thermo-refr active noise of the coating was taken into 
account as discussed in Sec. 14.8.2.61 

4.8.3.2 Sapphire substrate 

We have conducted the same kind of analysis for mirrors with Sapphire substrate. This 
study was made mainly for historical reasons as the idea of using sapphire substrates 
for the mirrors is now abandoned. In this case the dominant noise is the substrate 
thermo-elastic contribution. Some of the advantages of using Mesa bean have been 
already analyzed in [5D] for this particular noise source. Here we compute the various 
thermal noise contribution for finite size test mass and show the relative gain for each 
thermal noise employing a Mesa beam instead of a standard Gaussian. For 40 Kg 
sapphire substrates the minimum of thermal noise occurs for a mirror radius of about 



1 Bcnch is a simulation program available at http://emvogil-3.mit.edu/bench/ 
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4.8 Thermal noise: Gaussian VS mesa beam 



AdvLIGO Noise Curve with Gaussian and Mesa beam 




Frequency [Hz] 



Figure 4.6: Comparison between the Ad-LIGO expected sensitivity with Gaussian or 
mesa beams. 

16 cm. The corresponding mirror aspect ratio {2a/ H) is about 2.6. The total thermal 
noise reduction for Mesa beam is about a factor 1.55 around the minimum. 

4.8.4 Further consideration 
4.8.4.1 Finite size effect 

We first want to give a quantitative estimation of the finite test mass (FTM) effect 



on the different kinds of thermal noise. In Fig. 4.8 we plot the fractional correction 
provided by our analysis with respect to the formulas found in literature for Gaussian 
beams in the semi-infinite mirror case. For mirror radius of 18 cm, which provides the 
minimum of thermal noise, the correction is around 20% for substrate Brownian and 
less than 10% for the other contributions. 

Lovelace [98] performed the computation of simple scaling laws for thermal noise 
with arbitrary beam shapes and semi-infinite mirrors. He also analyzed the finite test 
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Figure 4.7: Ratios between the total thermal noise with Gaussian and mesa beams. 
The vertical lines correspond to the minima of the noise curves. 

mass correction provided by our results with respect to his derivation in the case of 
mesa beam and the conclusions are similar to the FTM effect for Gaussian beam. 
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Figure 4.8: FTM effects for Gaussian beam (1 ppm diff. loss). 
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4.8 Thermal noise: Gaussian VS mesa beam 



Noise 


GB 


MB 


LG50 


CB 
CT 


4.72 • 1(T 19 
1.99 • 1(T 19 


2.8 • 1(T 19 
1.17- 10~ 19 


3.2 • 1(T 19 
1.95 • 10" 19 



Table 4.2: Noise calculation for GW sensitivity at 100 Hz. Units cm/V 'Hz 



4.8.4.2 High order LG modes 

It has been recently proposed [69] the use of high order LG beams and excited modes 
of meas beams to reduce the thermal noise in interferometric GW detectors. The 
problem is that although these kind of beam provide a comparable reduction of the 
substrate thermal noises with respect to the mesa beam, they are sensibly worse than 
the mesa beam for the dominant coating thermal noises. The radial nodes of the higher 
order LG beam, cause rapidly variation of the elastic field near the coating surface. This 
produces an higher elastic energy density stored in the coating, which boosts the coating 
Brownian thermal noise, and higher thermal gradients, which are responsible for the 
thermoelastic noise. To substantiate these considerations we shown a comparison in 
the elastic energy distribution inside the mirror between a standard Gaussian beam, a 
mesa beam and a LG50 mode, all with 1 ppm of diffraction losses. 
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Figure 4.9: Elastic energy for a Gaussian beam with w = 6.85 cm. 



Because of the ripples in the elastic energy distribution, we expect the coating ther- 
moelatic noise for the higher order Laguerre Gauss beams to become large in the low 
GW signal frequency regime. The characteristic thermal diffusivity length rj- = \ 
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4.8 Thermal noise: Gaussian VS mesa beam 
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Figure 4.10: Elastic energy for a mesa beam with 6 = 11 cm. 
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Figure 4.11: Elastic energy for LG50 with w = 4 cm. 
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4.8 Thermal noise: Gaussian VS mesa beam 



increases as the sensitivity frequency decreases, in such a way that the temperature fluc- 



tuations involves a larger volume of the test mass. In Fig. 4.12 we show the comparison 
between the Gaussian, mesa, and LG50 beams with respect to the coating thermoelatic 
noise. At low frequency (~ 10 Hz) the LG50 noise is almost three times larger than for 
Gaussian beam, which remains a factor of two above the mesa beam noise. 
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Figure 4.12: Comparison between the coating thermo-elastic noise trends with fre- 
quency. 
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4.9 Optimized coating 



4.9 Optimized coating 

The limit sensitivity of interferometric gravitational wave antennas is set by the thermal 
noise generated in the dielectric mirror coatings. These coatings are currently made 
of alternating low-loss quarter- wavelength high/low index material layers. This design 
yields the maximum reflectivity for a fixed number of layers, but not the lowest noise 
for a prescribed reflectivity. 

Realization that much of the thermal noise in the coating comes from the material 
of the high index layer [95] has lead to investigation of optimizing the individual layer 
thicknesses while preserving the needed reflectivity. This motivated our recent investi- 
gation^] of optimal coating configurations, to guarantee the lowest thermal noise for a 
targeted reflectivity. This communication provides a compact overview of our results, 
involving truncated periodically-layered configurations^} Possible implications for the 
advanced Laser Interferometer Gravitational wave Observatory (LIGO) are discussed. 

All interferometers presently in operation use quarter-wavelength (QWL) designs, 
which are known (Bragg theorem) to be optimal, in the sense that they yield the largest 
reflectivity for any fixed number of layers (or, equivalently, the smallest number of layers 
for any prescribed reflectivity). However, QWL coatings do not yield the minimum TN 
for a prescribed reflectivity, and hence are not optimal for GW interferometers, where 
the quantity that should be maximized is the visibility distance of the instrument, 
which in its turn is limited by the coating thermal noise. 

This section addresses the above optimization problem, and is organized as follows. 
We outline our comprehensive working model for the multilayer mirror reflectivity and 
Thermal noise. Than we present an examples of optimized coatings, focusing on peri- 
odic stacked- doublet multilayers, and compare them with the standard QWL syntheses. 
Conclusions and recommendations follow in the end of this section. 

4.9.1 Multilayer Coating Reflectivity 

We consider a planar multilayer dielectric coating composed of alternating homogeneous 
layers (with variable thickness) of silica (S1O2) and tantala (To^Os), illuminated by a 

1 In collaboration with the Waves Group, Department of Engineering, University of Sannio, 
Benevento, Italy. 

2 More general analysis involving nonperiodic genetically-engineered coatings has been pre- 
sented in [55] 
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normally-incident plane wave, with implied time-harmonic dependence. Referring to 
Fig. |4. 13 we consider M layers, M + 1 interfaces, and four dielectric media, including 
the left and right semi-infinite media a and b (vacuum and substrate). 
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Figure 4.13: Multilayer dielectric structure. 

The incident and reflected fields are considered at the left of each interface. The 
overall reflection response, Ti = E1—/E1+, can be obtained recursively in a variety 
of ways, such as the propagation matrices, the propagation of the impedances at the 
interfaces, or the propagation of the reflection responses (IQOj . The elementary reflec- 
tion coefficients pi from the left of each interface are defined in terms of the refractive 
indices as follows: 



Pi 



Tlj-l - Tij 



J 



1, 



M + l 



(4. 



rij-i + no- 
where pi = (ra - ni)/(n a + m) and p M +i = {n M - n h )/(n M + n b ). The for- 
ward/backward fields at the left of interface i are related to those at the left of interface 
j + 1 by: 
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where kjlj is the phase thickness of the j — th layer, which can be expressed in 
terms of its optical thickness njlj and the operating free-space wavelength by kjlj = 
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2n(rijlj)/X. The reflection responses Tj = Ej_/Ej + will satisfy the recursions: 



-47T JJ 



r _ Pi + r i+ie —A j = M, M — 1, . . . , 1 (4.89) 

i + Pi r J+ie - 47r — 

and initialized by Tm+i = pm+i- The mirror reflectivity is given by V = Y\. It is useful 
for the following analysis to introduce the scaled optical length Zj = rijlj/X. 

The above formalism can be easily generalized |100] to the case of oblique incidence, 
via the introduction of (polarization dependent) "transverse" wave-impedances and 
refractive indexes. 

4.9.2 Multilayer Coating Thermal Noise 

Since the function we want to minimize is the coating thermal noise, we need the most 



accurate description of the thermal noise for a multilayer structure. Eq. (4.49) has 
been generalized by Harry [101J in order to incorporate coating anisotropy into the 
thermal noise. A multi layer structure composed of alternating layers of homogeneous 
materials exhibits a mechanical behavior proper of a transversely isotropic material in 



which the independent elastic coefficient in the stiffness matrix (4.19) are reduced to 
five by the rotational symmetry around the axis orthogonal to the layers. Ki and Y± 
are the Young modulus for stress causing strains entirely within the plane parallel to 
the coating layers or perpendicular to this one. There are also two Poisson ratios, uu for 
stresses and strains both with the plane parallel to the coating layers, and is± for when 
either the stress or the strain is perpendicular to the coating layers. The displacement 
induced noise is given by 

\JlTUJ WI S JJ 

where the loss angle 9^// depends on the coating thickness d, on the beam radius w, 
and on the loss angles in the parallel and perpendicular direction with some coefficients 
dependent on the elastic coefficient of the coating and substrate 
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.coat f^±_ 2 ^ Y \\ 

e ff . a;,,, v, \\ i _ ,,2 



ttwYjI U-i/f Y±(l -v 2 s ){l -v\\) 
Y»v 1 _{l-2v s ) 
+ (1-^(1-^)^-^ 



+ w-ifjd-,.) 011 J (491) 

In the paper |101j . Y±, Y\\ , i/m , <fi±, 4>\\ , are calculated from the values of the 
isotropic materials that makes up the layers of the coating in the following way 



Y ± = (d 1 + d 2 )/(d 1 /Y 1 + d 2 /Y 2 ) 

Y {] = (Y 1 d 1 + Y 2 d 2 )/(d 1 + d 2 ) 

<p x = y±(0idi/yi + 2 d 2 /r 2 )/(rf 1 + d 2 ) 

0,1 = (Fi^irfi + Y 2 <j> 2 d 2 )/[Y\\ (di + da)] 

u ± = {v 1 Y 1 d l + v 2 Y 2 d 2 )/{Y 1 d l +Y 2 d 2 ) 

v\\ = F(Y ll Y 2l vi,v 2 ,d ll d 2 ) 



v» solution of 



(l + ^i)(l-2^) 



^ 2 F 2 d 2 
(l + ^ 2 )(l-2^ 2 ) 



>ll + l)(2^i| | _(i_ I/|| )y ± ) 



When we first looked at the problem of coating optimization, we were not completely 
satisfied with the derivation of the averaged mechanical properties of the multi-layer 
coating quoted above. Therefore we put some effort in the rigorous analysis of a multi- 
layer structure and the results are presented in the following section. 

4.10 Averaged elastic coefficients for multilayered optical 
coating 

Heterogeneity of a material or structure can be caused by two main reasons. 

1. Non-uniformity of certain physical characteristic (density, elastic modulus, con- 
ductivity, etc. ). Two or multi-phase composites are typical examples of this type 
of material behavior. 
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2. An additional source of heterogeneity is a geometrical one. A multi-layer structure 
like the dielectric optical coating used in gravitational waves interferometers is an 
example of this effect. 

The replacement of heterogeneous media by an homogeneous continuum, which is 
characterized by certain effective constitutive equations, is the basic instrument for the 
effective media theory (EMT), which is a well developed subject in composite mechanics 
[102 . The purpose of this section is to give a mathematical introduction to the prob- 
lem of the calculation of mechanical properties of heterogeneous solids with periodic 
microstructure. The averaging procedure for mechanical fields and the definitions for 
effective constitutive parameters (effective elastic moduli, for instance) are described 
in the following for a bilayered composite medium. We provide a rigorous derivation 
of the coating effective mechanical parameters which is very important in the study of 
coating optimization for GW interferometers. 

Let us consider an elastic solid which consists of a periodic array of two layers 
with thicknesses d\ and cfe, respectively in welded contact and stacked in the X3 = z 
direction. It means that the components, Cijki, of the elastic moduli tensor, which link 
strains e^i, and stresses, Cjj are periodic functions, Cijki(z + d\ + d%) = Cijkl{z)- Let 
us define the average over the unit cell, composed dy a bilayer structure, as 

(X) = dl * ± d f = 5l X + 5 2 X 
«1 + 0L2 

One can introduce the definition of the effective moduli tensor C*- kl as the coeffi- 
cients which link the components of effective fields variables defined by averaging the 
strain and stress tensor and taking into account the continuity conditions at the inter- 
face between the two layers. We skip all the mathematical technicalities, referring the 
reader to |102| 1103] and we report only the results of the homogenization process 
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It is interesting to point to the direct averaging law for the composite planar shear 
modulus C{ 2 i2 and the inverse averaging law for the transverse ones: Ci3i3>C2323 - 
The constitutive relations for each of the two assumed isotropic materials are given 



by (4.23) as functions of the two Lame coefficients. In a cartesian coordinate system, 
the components of the stiffness tensor are 



(4.93) 



Let us remind in Tab. 4. 10 the relation between the Lame coefficients and the Young 
and Poisson coefficient for an isotropic material 



Y 



V ~ 2{l+v) 



^t(3A+2^)~ 
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Yv 

(l+v){l-2v) 



Y 



A 

2(\+ti) 



Table 4.3: Elastic coefficients 

Since both materials are assumed isotropic, the resulting layered medium is expected 
to posses transverse isotropy with isotropy confined in the x\ — %2 plane. Even in 
an anisotropic medium, the axial components of the stiffness tensor are frequently 
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characterized by the parameters of uniaxial tension, namely Young modulus Yj, and 
the Poisson ratio vji, for the tension in the Xi direction: 



Yi = — and Vji = — — ,ajj = 0,j i (4.94) 



These parameters can be easily calculated for the homogenized compound material 
by means of the components of the effective stiffness tensor C*- kl 



* iijj kkkk iikk jjkk , . „_\ 

^jjjj^kkkk y^jjkk) 

*i ^iiii iijj ji ^iikk^kii 

i^j,i^k,j^k (4.96) 

Planar, Y* = Y 2 * = Y\\ an d transverse Y 3 * = Y±, Young moduli, share moduli and 
Poisson ratios of binary multilayered coating composed of S1O2 and TCI2O5 are shown 
in fig. REF as functions of the Ta20^ layer thickness fraction. The explicit expressions 
of the elastic coefficient of the multilayer structure as function of the elastic coefficients 
of the two material are 



Y 1 {l-vl)S 1 + Y 2 {l-vl)S 2 

Y 1 Y 2 [Y 1 (l-u 2 )S 1 + Y 2 (l-iy 1 )5 2 ] 



v * _ v * _ 1 V 2/1 1 J 2 \ ^ v\v 2 )u\u 2 T i 2 \i- - vpu2 ( . Q7 N 



Y * Y 2 2 (l -v x - 2v 1 2 )5 1 5 2 + F 2 2 (l - v 2 - 2v%)5 1 5 2 + Y X Y 2 [(1 - v 2 )8\ + 4iWi& + (1 - Vl )b\\ 

* _ Y lVl {l- v l)5 1 + Y 2 u 2 {l-ul)5 2 

^(1-^ + ^(1-^ (4 - 98) 

* = Yjj [(l - ^2)^1^1 + (l - ^1)^2] 

Yftl - v 2 - 2vl)b x b 2 + F 2 2 (l - v x - 2^)6^2 + Y X Y 2 [(1 - v 2 )5\ + Av x v 2 5 x 5 2 + (1 - i^] 
YY 2 



2[Y 2 (1 + j/ 1 )£ 1 +Y 2 (1 + ia,)<5 2 ] 



^ _ Yifr Y 2 S 2 



2(l + i^) 2(1 + iy 2 ) 



By convention, the 5 elastic constants in transverse isotropic constitutive equations 
are the Young's modulus and Poisson ratio in the x-y symmetry plane, Yy and z/y, 
the Young's modulus and Poisson ratio in the z-direction, Yj_ and isj_, and the shear 
modulus in the z-direction G±. The compliance matrix takes the form 
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Figure 4.14: Averaged elastic coefficients using material parameters of App. A. 3 
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Let us now expand the exact expressions of the Young and Poisson averaged coeffi- 
in power series of the two material's Poisson coefficient v\ and 



cient, Eqs. fl4.97[ ),(|4. 



Y {{ = YtSi + Y 2 5 2 + 0{u 2 ) 

Y ± = ™ 

Y 2 5i + Y^ 2 V ; 



(4.99) 



' y ll 



Y\v\b\ + Y 2 v 2 5 2 



+ 0(v 2 



Y X S X + Y 2 5 2 
YxY^xdx + YxY 2 u 2 S 2 



+ 0(v 2 



{Y 1 6 1 + Y 2 5 2 )(Y 2 6x + Yx5 2 ) 

These more manageable expression^] offer excellent approximation for typical LIGO 
mirrors coating: for the standard quarter wavelength design, the difference between the 
approximated values and the exact values is within 0.5%. 



4.11 Thermal noise for an anisotropic coating 

The thermal noise for a semi-infinite mirror with an anisotropic coating is given by |^] 



Sx(f) 



2K B Tl-u 2 coat 



3 , 

7T2/ 



wY x 



J eff 



(4.100) 



lr The firsts two expression reproduce the results of |101j . whereas the other two are different. 
2 This is a corrected form of the equation (8) in |101j . 
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4.11 Thermal noise for an anisotropic coating 



h coat d 1 ( / Y s _ 2v\Y s Y\\ 
eff Y ± Ul - i/2 Y ± {1 - z/|)(l - i/n) 



Yi|i/j_(l - 2i/ s ) 

' r '-f±) 



(1-^(1-^) ^11 

Yiiyi(i + t/ a )(i-2t/ 8 ) 2 
y,(i-^)(i-i/ a ) 



(4.101) 



In the limit of isotropic coating this formula reduces to (4.49 1 . Since the loss angle 
is introduced in the theory of elasticity as an imaginary part of the elastic coefficient of 
the material, it is natural to obtain the values of <f>?, and 4>^_ from the isotropic values 
01) 02) using the same averaging rules as the Young moduli. Therefore we have 

5 =^ and 01 =<^ (4 , 02) 
and the expansion in the Poisson coefficients gives 



^11 ~ VA , vx + °K V aIld 4>± = VAX _j_ \r i f + (4.103 

11 YiO\ + F2O2 *l<Pl<32 + ^20201 

We now have all the ingredient to express the coating thermal noise as function 
of the characteristic elastic coefficient of the substrate and coating materials, which 
for the present analysis will be treated as constants, and the thickness, 81,82, and the 
number of doublet, of the layers. The function Sx(8i, 82, N^) results in a complicated 
expression which is not particulary illuminating and its analytical form will be omitted 
in this presentation. However we can give a very good approximation of this function 
by a first order Taylor expansion around the point which correspond to the standard 
design, z\ = z 2 = 1/4, where z is the scaled optical path z = n8/X. 

Sx (f) = 2K BTl-^ s Na (221Q -lo Zi + 1A5W -n Z2) (4 1Q4) 

7T2/ WY S W 



Fig. 4. 15 left panel) shows the accuracy level of this linear approximation. If we are 



looking for an analytical approximated expression we can take the limit of vanishing 



151 



4.12 Coating thickness optimization 



Poisson ratios in the general expression of the noise and then take the linear terms in 
the two optical thickness Z\ , z 2 



Sx(f) 



1K B T\-vl N d X 



tt! / wY s \m \Y S Yj Zl + n 2 \Y S + Y 2 1 1 



(4.105) 



From Fig, 4.15 right panel) we can easily see that even if the expression above is 
very manageable, it provides a rough approximation of the exact formula and it leads 
to error above 30% in the interesting region. 





Figure 4.15: Left: Fractional error using expression 4.104 Right: Fractional error for 



the expression 4.105 



4.12 Coating thickness optimization 

Our collaborators from the University of Sannio, Benevento, developed an analysis 
based on genetic algorithm (GA) to address the problem of finding the best coating 
design for minimizing the thermal noise. We briefly review their conclusions, in order 
to validate our simple direct optimization. Genetic algorithms (GAs) 

GAs have proven as effective tools for synthesizing general (e.g., multi-dielectric) re- 
flective coatings, featuring several heterogeneous (e.g., technological) design constraints, 
with multiobjective (i.e., reflectivity and thermal noise) optimization targets. The gen- 
eral analysis led to the result that, as the number of generations is increased, the 
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4.12 Coating thickness optimization 



centermost part of the GA-optimized coatings exhibits a neat tendency toward a cas- 
cade of almost identical doublets, whose (total) optical lengths cluster around A/2. The 
above findings led us to investigate the performance of coatings composed of cascaded 
identical multiplets with the constraint z\ + zi = 1/2. This constraint reduces the 
number of independent variables to two, z\ and Nd, the optical thickness of the high 
index material and the number of doublets. 



Fig. 4.16 shows the noise function Sx(zi,Nd) scaled to the Ad-LIGO baseline value, 



and 1 — |r| (zi, Nd) in the range of interest. 
We followed the strategy: 

• Calculate the coating transmission^] for the baseline design, z\ = Z2 = 1/4, Nd 



19, using Eq.(|4.89l): 1 - \T\ 2 = 8.3ppm 



In the plane (z±,Nd), draw the curve of equal transmission constrained by the 
above value. 

Plot the noise function Sx(z%, Nd) in the same graph and find the point on the 
above curve which correspond to the minimum noise. This correspond to the 
optimized design (zi, Nd) op t- 



All the procedure is synthesized in Fig. 4. 17 The optimized design (for the known 



values of material losses, elastic moduli, etc.) is given by Nd = 23, z\ = 0.153 and Z2 = 
0.347; for the same optical transmission of the A/4-design, it provides a 20% of reduction 
in noise spectral density. Under the simplest assumption where the GW sources are 
distributed homogeneously /isotropically throughout space, without interferometer re- 
optimization, this may boost the event rate by some 30%. 

The optimization of the coating structure to minimize the thermal noise is therefore 
a very promising idea, which in principle could stand beside the mesa beam proposal, 
increasing further the detector sensitivity. Coating thickness optimization should be 
considered as almost mandatory to minimize coating noise, yielding in all cases^ja sub- 
stantial increase(> 30% at 100 Hz) in the expected event rate, as compared to the 

1 In the absence of (optical) losses. 

2 In this thesis we have focused the attention on the Brownian coating thermal noise, but 
on-going researches [57] have analyzed the total coating noise budget and provided the corre- 
sponding optimized design. 
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(a) Thermal noise contour plot scaled to the Ad-LIGO baseline de- 
sign. 
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(b) Iso-transmission curves. 

Figure 4.16: The constraint z\ + Z2 = 1/2 allows a simple optimization by inspection 
of these plots . 
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zi 



Figure 4.17: The blue curve corresponds to a transmission of 8.3ppm, the green point 
corresponds to the Ad-LIGO coating standard design, and the red point corresponds 
to the optimized design. 

QWL design. A prototype of optimized coating has been designed by our collaboration 
and is scheduled for testing at Caltech in the TNI (Thermal Noise Interferometer) fa- 
cility. LMA is responsible for manufacturing the sample, using the ion beam deposition 
techniques. 
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Chapter 5 

Conclusions 



This thesis presents a collection of different researches on non-standard optics in view 
of enhancing the performances of the Advanced Gravitational waves interferometric 
detectors. 

Although primarily based on simulations, calculations and theoretical arguments, 
this work has been in continuous feedback with the experimental side of the thermal 
noise minimization effort. This work has provided the theoretical understanding and 
guidance for the experimental side, while getting back continuous and progressive vali- 
dation of the predictions and prescriptions. The end result was a very good example of 
synergy, and mutual enrichment, for me and my introduction to laboratory work, and 
for my colleagues that worked on my prescriptions and numbers. 

In the design of the next generation Gravitational Waves interferometric detectors, 
mirror's thermal noise plays a crucial role. The mirrors coating thermal noise is ex- 
pected to be the limiting factor in the GW signal frequency region of highest sensitivity. 

The quest for increasing the event rate in the observational band has motivated 
the exploration of various techniques for reducing the mirror thermal noise toward 
the development of second-generation detectors, such as Advanced-LIGO [16]. With 
specific reference to the coating Brownian noise (dominant in the current baseline design 
featuring fused-silica test masses), use of improved (low- mechanical-loss) materials, 
geometric optimization of the coating design, and flat-top (commonly referred to as 
"mesa") beams seem the most promising. The latter option, intuitively motivated 
by the potential capability of a mesa beam (MB) of better averaging the thermally 
induced mirror surface fluctuations as compared to a standard Gaussian beam (GB), 
has been numerically proved in this thesis to yield significant reductions in the overall 
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thermal noise. For the first time we calculated the coating types of thermal noise 
reduction in using mesa beam and enlightened the issue of mirror geometric aspect 
ratio optimization. We demonstrated that the coating induced displacement noise can 
be reduced by a factor of ~ 1.8 which increases the estimated binary neutron star 
inspiral range Q seen by Ad-LIGO, from 175 Mpc to 225 Mpc. 

During this work we implemented the different kinds of thermal noise calculations 
for finite size test mass and arbitrary cylindrically-shaped laser beams, in a simple 
Mathematica® notebook called Thermal Noise Notebook (TNN). This package is avail- 
able as a very simple tool for the estimation of the thermal noise contributions taking 
into account both the beam and the mirror geometry. We illustrated the importance 
of uniform sampling of the mirror surface to reduce thermal noise and the limitation 
brought by the use of excited modes with nodes on the mirror surface, which cause 
non-negligible gradients in the elastic and thermal fields in the test mass and worsen 
the thermal noise performance. 

In Sec ]3.4| we further developed the theory of mesa beam and derived some unex- 
plored analytical results for the dual configuration, flat mesa beam and concentric mesa 
beam, which are of interest in view of the implementation of this non-standard optics 
in real systems. The expression of the beam width, divergence and beam propagation 
parameter, M 2 , are derived in analytical form as functions of the characteristic param- 
eters of the mesa beams, thus allowing a straightforward application of the well known 
ABCD formalism for the propagation of optical beams trough paraxial systems. 

We also analytically proved a new duality relation between optical cavities with non- 
spherical mirrors. This derivation provides a unique mapping between the eigenvalues 
and eigenvectors of two cavities whose mirrors shapes are related by a simple relation. 
This duality allows the direct application of beam property calculations performed in 
a case to geometries of the other configuration. 

The growing interest in the GW community in the mesa beam idea, has led to 
the development, by our group, of a Fabry-Perot cavity prototype with non-spherical 

1 Thc effective range conventionally measures the sensitivity of the interferometers to signals 
arising from binary systems composed of 1.4 solar mass neutron star pairs having a signal-to- 
noise ratio of 8 : 1 or better, and is computed by averaging over both the antenna sensitivity 
pattern and over all possible orientations of the binary systems. 
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"Mexican hat" (MH) profile mirrors. One of the main task in this thesis has been 
the development and testing of optical simulation packages based on FFT routines 
or on Huygens paraxial integral approach. The FFT based program can simulate a 
Fabry-Perot optical cavity that includes non spherical mirrors and / or "realistic" mirror 
deformations. Many different cavity parameters or imperfections can be modelled as 
well as mirror misalignment and input beam coupling sensitivity. These programs have 
been successfully used to provide a theoretical framework for the experimental work, 
i.e. the fundamental and higher modes sensitivity to cavity misalignment and mirror 
imperfections, for our Mexican hat mirror cavity, for the large scale Ad-LIGO FP 
cavities or any interferometer to be designed. Some examples of their application are 
given in this thesis as well as a discussion on the accuracy and limitations they provide. 

The research on the optical coatings for next generation detectors requires the most 
possible accurate model of the coating elastic and thermal parameters. In this thesis 
we present a rigorous, never used in this problem (to our best knowledge), derivation 
of the coating multi-layer effective elastic coefficients as function of the constituents 
individual properties. Some simplified expression are also derived, as leading term of 
the expansion of the exact formulas, which are in agreement with published results. This 
analysis is of fundamental importance when dealing with the comparison of different 
coating design because the elastic coefficients of the coating enter the thermal noise 
evaluation formulas. 

The last part of the thesis has been dedicated to the study of the coating thermal 
noise reduction by a modification of the geometric structure of the multilayered coating. 
The scope of the project has been to find an alternative coating design respect to 
the standard quarter-wavelength one, which, for a prescribed transmittance, gives the 
minimum thermal noise. We found that the optimized design decreases the mirror 
thermal noise of about 10%. Under the simplest assumption where the GW sources 
are distributed homogeneously /isotropically throughout space, without interferometer 
re-optimization, this may boost the event rate by some 30%. A prototype of optimized 
coating has been designed by our collaboration and is scheduled for testing at Caltech 
in the TNI (Thermal Noise Interferometer) facility. 

All these research topics are currently very actively investigated by our collabora- 
tions |68j . The end result of our studies will likely be the extension of the reach of GW 
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observatories by a factor of a few, and of the event rate by a factor of several, maybe 
more than an order of magnitude. 
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Appendix 



A.l Physical components of tensors in orthogonal coordi- 
nate systems 

The physical components of a tensor are given by the components of the tensor with 
respect to unit basis vectors. In an orthogonal coordinate system, the metric tensor is 
diagonal and is conventionally written as 



For example, the physical component of any vector in a given direction is merely 
the projection of that vector onto that direction, i.e., the inner product of the vector 
with a unit vector in the given direction. This is simply the component of the vector 
in a basis having a unit basis vector pointing in that direction. 

The unit contravariant vectors tangent to the three coordinate lines are 




(A.2) 



The index enclosed in parentheses means a non tensorial index. In terms of the co- 
variant/contravariant and metric components, the physical components of an arbitrary 
vector ,in an orthogonal coordinate system, are given by 



= hf , g u = ~2 , mot summed, and gij = = i / j 



(A.l) 




no summation over indices 



(A.3) 



and for a mixed tensor 



T(i..jk..l) 




no summation over indices 



(A.4) 
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A.l Physical components of tensors in orthogonal coordinate systems 



The covariant derivative of a scalar function which reduces to the partial derivative, 
is a covariant vector the gradient of <£. The physical components of the gradient 
of a scalar function are given by 



grad&(i) 



1 d<S> 

hj dx % 



(A.5) 



The covariant derivative with respect to x l of a covariant vector A 1 summed on % is the 
divergence of the vector A and can be written in terms of the physical components of 
the vector, A(i) 



divA 



1 d 



(g 1/2 A* 



1 d 



1/2 



A(i) 



(A.6) 



g 1 / 2 dx' 1 g 1 / 2 dx l 

where g denotes the determinant of the metric and in orthogonal coordinates is 
expressed as g 1 ^ 2 = /ti 712/1-3. 

The Christoffel symbols in orthogonal coordinates are 



r 0, 



r i = V -n ir ( ^± + _ dffjfc \ 

jk Z^ 2 9 \dx k dxi dx r J 



1 dh t 

hi dxi ' 



when i,j, k are all different; 

when i = j = k, 

i = j / k,i = k ^ j; 



The physical components of the covariant derivative of a vector field are 



A(i;j) 



-A: 



_d_ 

dxi 



(A.7) 



(Ai 



Using equation (A. 8) and (A.7) the physical components of the strain tensor can 



be written as function of the physical components of the displacement vector 



e(ii) 



d ( u(i) 
dx % \ hi 



1 u(m) d . 
+ TrT~n / — ^ — [h i ) no summation on 1 



2h 2 ^ h m dx r 



hi d f u{i) \ hj d f u(j) 



hj dxi \ hi 



+ 



hj dx 1 



(A.9) 

no summation on i,j, i ^ j (A. 10) 
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A. 2 Equilibrium equations 



The expansion in terms of the physical components of the displacement vector 



can be written using (A. 6) 



6 



1 



hih 2 h 3 



d(h 2 h 3 u(l)) d(h 1 h 3 u(2)) dihxhiui?,)) 



dx 1 



dx 2 



dx 3 



(All) 



A. 2 Equilibrium equations 

The conditions for static equilibrium in a linear elastic material are determined from 
the force balance equation in a generic coordinate system 



(A.12) 



where a™ are the stress tensor components,/' are the external body forces per unit 



mass and p is the density of the material. Using the Eq.(A.8) and (A. 7) we can express 
the equilibrium equation in terms of physical components 



E 



hi d 



V9 



a(ij)j+^^-r—j—T t j k a(Jk) + pf(i) = no summation on i (A. 13) 



0,k 



Using the constitutive equations for a linear isotropic material Eq.(4.24), and the 



definition of the strain tensor in terms of the displacement vector Eq.(4.16) we have 



that the equilibrium equations Eq. (A.12) becomes 



(A + fl)Uk;ki + P u i:jj + Pfi = 



(A.14) 



In an orthogonal coordinate system, the components of V(V • u) can be expressed 
in terms of physical components by the relation 



[V(V-u)], 



LA. 

hi dx' 1 



E 



1 d 



hk 



u(Jfe) 



(A.15) 



; ^/gdx k 

and the components of the second term becomes quite complicated in terms of the 
physical components of the displacement vector 
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d 2 {h iU (i)) 




- 2 £r 



d(h m u(m)) 



,d (hju{i)) 
dx m 



\ *P 33 JP *j 




(A.16) 
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A. 3 Mirror's elastic and thermal coefficients 



A. 3 Mirror's elastic and thermal coefficients 



Parameters 


Fused Silica 


Sapphire 


Ta 2 5 


Density p (g/cm s ) 


2.2 


4 


6.85 


Young modulus Y (erg /cm 3 ) 


7.2 • 10 11 


4 • 10 12 


1.4- 10 12 


Poisson ratio v 


0.17 


0.29 


0.23 


Loss angle <p 


5 • 10~ 9 


3 • 10~ 9 


4.5 • 10~ 4 


Linear thermal expansion a (K^ 1 ) 


5.5 • 10~ 7 


5 • 10~ 6 


3.6- 10~ 6 


Specific heat per unit mass C (erg/(g K)) 


6.7- 10 6 


7.9- 10 6 


3.06 • 10 6 


Thermal conductivity K (erg / (cm s K)) 


1.4 • 10 5 


4- 10 6 


1.4 • 10 5 


Refraction index (at A = 1.064/im) n 


1.46 


1.75 


2.06 



Table A.l: These parameters are referred to a temperature of 300 K. In reference [92] 
the thermal conductivity of Ta^O^ is chosen equal to that one of crystalline sapphire, 
but we think it would be more appropriate to take a value more similar to that of 
amorphous fused silica, being a ion-beam sputtered thin layer [104J. Moreover the 
thermal conductivity of most thin films is found [105] to be orders of magnitude lower 
than that of the material in bulk form. How this affects the thermal noise evaluation 
for the GW interferometers is under investigation. 
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